# multi_physics/piezo_elastodynamic.py¶

Description

The linear elastodynamics of a piezoelectric body loaded by a given base motion.

The generated voltage between the bottom and top surface electrodes is recorded and plotted. The scalar potential on the top surface electrode is modeled using a constant L^2 field. The Nitsche’s method is used to weakly apply the (unknown) potential Dirichlet boundary condition on the top surface.

Find the displacements , the potential and the constant potential on the top electrode such that:

where is the matrix of elastic properties under constant electric field intensity, the piezoelectric modulus, the permittivity under constant deformation, a penalty parameter and the external circuit resistance (e.g. of an oscilloscope used to measure the voltage between the electrodes).

## Usage Examples¶

Run with the default settings, results stored in output/piezo-ed/:

sfepy-run sfepy/examples/multi_physics/piezo_elastodynamic.py


The define() arguments, see below, can be set using the -d option:

sfepy-run sfepy/examples/multi_physics/piezo_elastodynamic.py -d "order=2, ct1=2.5"


View the resulting potential on a deformed mesh (2000x magnified):

sfepy-view output/piezo-ed/user_block.h5 -f p:wu:f2000:p0 1:vw:wu:f2000:p0 --color-map=inferno


source code

r"""
The linear elastodynamics of a piezoelectric body loaded by a given base
motion.

The generated voltage between the bottom and top surface electrodes is recorded
and plotted. The scalar potential on the top surface electrode is modeled using
a constant L^2 field. The Nitsche's method is used to weakly apply the
(unknown) potential Dirichlet boundary condition on the top surface.

Find the displacements :math:\ul{u}(t), the potential :math:p(t) and the
constant potential on the top electrode :math:\bar p(t) such that:

.. math::
\int_\Omega \rho\ \ul{v} \cdot \ul{\ddot u}
+ \int_\Omega C_{ijkl}\ \veps_{ij}(\ul{v}) \veps_{kl}(\ul{u})
- \int_\Omega e_{kij}\ \veps_{ij}(\ul{v}) \nabla_k p
= 0

\int_\Omega e_{kij}\ \veps_{ij}(\ul{u}) \nabla_k q
+ \int_\Omega \kappa_{ij} \nabla_i \psi \nabla_j p
- \int_{\Gamma_{top}} \kappa_{ij} \nabla_j p n_i q
+ \int_{\Gamma_{top}} \kappa_{ij} \nabla_j q n_i (p - \bar p)
+ \int_{\Gamma_{top}} k q (p - \bar p)
= 0

\int_{\Gamma_{top}} \kappa_{ij} \nabla_j \dot{p} n_i + \bar p / R = 0 \;,

where :math:C_{ijkl} is the matrix of elastic properties under constant
electric field intensity, :math:e_{kij} the piezoelectric modulus,
:math:\kappa_{ij} the permittivity under constant deformation, :math:k a
penalty parameter and :math:R the external circuit resistance (e.g. of an
oscilloscope used to measure the voltage between the electrodes).

Usage Examples
--------------

Run with the default settings, results stored in output/piezo-ed/::

sfepy-run sfepy/examples/multi_physics/piezo_elastodynamic.py

The :func:define() arguments, see below, can be set using the -d option::

sfepy-run sfepy/examples/multi_physics/piezo_elastodynamic.py -d "order=2, ct1=2.5"

View the resulting potential :math:p on a deformed mesh (2000x magnified)::

sfepy-view output/piezo-ed/user_block.h5 -f p:wu:f2000:p0 1:vw:wu:f2000:p0 --color-map=inferno
"""
from functools import partial

import numpy as nm

from sfepy.base.base import output
from sfepy.discrete.fem.meshio import UserMeshIO
from sfepy.mesh.mesh_generators import gen_block_mesh
from sfepy.homogenization.utils import define_box_regions

def post_process(out, problem, state, extend=False, pcs=None):
"""
Calculate and output strain, stress and electric field vector for the given
displacements and potential.
"""
from sfepy.base.base import Struct

ev = problem.evaluate
strain = ev('ev_cauchy_strain.i.Omega(u)', mode='el_avg', verbose=False)
stress = ev('ev_cauchy_stress.i.Omega(m.C, u)', mode='el_avg',
copy_materials=False, verbose=False)

out['cauchy_strain'] = Struct(name='output_data', mode='cell',
data=strain)
out['cauchy_stress'] = Struct(name='output_data', mode='cell',
data=stress)
out['E'] = Struct(name='output_data', mode='cell', data=E)

top = problem.domain.regions['Top']
p_top = state['p'].get_state_in_region(top)
# = state['pc'](), but we want to test .get_state_in_region()
pc_top = state['pc'].get_state_in_region(top)

output('pc:', pc_top)
output('|p - pc|_top:', nm.linalg.norm(p_top - pc_top))
if pcs is not None:
pcs.append(pc_top[0, 0])

return out

def plot_voltage(problem, state, pcs=None):
import os.path as op
import matplotlib.pyplot as plt

ts = problem.get_timestepper()

fig, ax = plt.subplots()
ax.plot(ts.times, pcs)
ax.set_xlabel('$t$ [s]')
ax.set_ylabel(r'$\bar p$ [V]')

fig.tight_layout()
fig.savefig(op.join(problem.output_dir, 'voltage.pdf'))

def define(
dims=(1e-2, 1e-2, 5e-3),
shape=(5, 11, 21),
order=1,
amplitude=0.0000001,
ct1=1.5,
dt=None,
tss_name='tsn',
tsc_name='tscedl',
ls_name='lsd',
active_only=False,
save_times='all',
output_dir='output/piezo-ed',
):
"""
Parameters
----------
dims: physical dimensions of the block mesh
shape: numbers of mesh vertices along each axis
order: the FE approximation order
ct1: final time in min(dims) / "longitudinal wave speed" units
dt: time step (None means automatic)
tss_name: time stepping solver name (see "solvers" section)
tsc_name: time step controller name (see "solvers" section)
ls_name: linear system solver name (see "solvers" section)
save_times: number of solutions to save
output_dir: output directory
"""
dim = len(dims)
assert dim == 3

# A PZT 5-H material, Voigt notation, strain - electric displacement form.
epsT = nm.array([[1700., 0, 0],
[0, 1700., 0],
[0, 0, 1450.0]])
dv = 1e-12 * nm.array([[0, 0, 0, 0, 741., 0],
[0, 0, 0, 741, 0, 0],
[-274., -274., 593., 0, 0, 0]]) # C / N = m / V

# Convert to stress - electric displacement form.
CEv = nm.array([[1.27e+011, 8.02e+010, 8.47e+010, 0, 0, 0],
[8.02e+010, 1.27e+011, 8.47e+010, 0, 0, 0],
[8.47e+010, 8.47e+010, 1.17e+011, 0, 0, 0],
[0, 0, 0, 2.34e+010, 0, 0],
[0, 0, 0, 0, 2.30e+010, 0],
[0, 0, 0, 0, 0, 2.35e+010]])
ev = dv @ CEv
epsS = epsT - dv @ ev.T

# SfePy: 11 22 33 12 13 23
# Voigt: 11 22 33 23 13 12
ii = [0, 1, 2, 5, 4, 3]
ix, iy = nm.meshgrid(ii, ii, sparse=True)
CE = CEv[ix, iy]
e = ev[:, ii]

eps0 = 8.8541878128e-12
kappa = epsS * eps0

# Longitudinal and shear wave propagation speeds.
mu = CE[-1, -1]
lam = CE[0, 0] - 2 * mu
rho = 7800
cl = nm.sqrt((lam + 2.0 * mu) / rho)
cs = nm.sqrt(mu / rho)

# Element size.
L = nm.min(dims)
H = L / (nm.max(shape) - 1)

# Time-stepping parameters.
if dt is None:
# For implicit schemes, dt based on the Courant number C0 = dt * cl / H
# equal to 1.
dt = H / cl # C0 = 1

t1 = ct1 * L / cl

# Time history record of pc.
pcs = []
_post_process = partial(post_process, pcs=pcs)
_plot_voltage = partial(plot_voltage, pcs=pcs)

def mesh_hook(mesh, mode):
"""
Generate the block mesh.
"""
mesh = gen_block_mesh(dims, shape, 0.5 * nm.array(dims),
name='user_block', verbose=False)
return mesh

elif mode == 'write':
pass

filename_mesh = UserMeshIO(mesh_hook)

bbox = [[0] * dim, dims]
regions = define_box_regions(dim, bbox[0], bbox[1], 1e-5)
regions.update({
'Omega' : 'all',
})

materials = {
'inclusion' : (None, 'get_inclusion_pars')
}

fields = {
'displacement' : ('real', 'vector', 'Omega', order),
'potential' : ('real', 'scalar', 'Omega', order),
'constant' : ('real', 'scalar', 'Top', 0, 'L2', 'constant'),
}

variables = {
'u' : ('unknown field', 'displacement', 0),
'v' : ('test field', 'displacement', 'u'),
'p' : ('unknown field', 'potential', 1, 1),
'q' : ('test field', 'potential', 'p'),
'pc' : ('unknown field', 'constant', 2, 1),
'qc' : ('test field', 'constant', 'pc'),
}

materials = {
'm' : ({
'C': CE,
'e' : e,
'kappa' : kappa,
'rho': rho,
'penalty': 1,
'iR' : 1.0 / (15e6 * dims[0] * dims[1]), # 1 / (R * top_area).
},),
}

integrals = {
'i' : 2 * order,
}

def get_ebcs(ts, coors, mode='u'):
y = coors[:, 1]
k = 2 * nm.pi / dims[1]
shift = nm.pi / 3
omega = cl * k
time = ts.time
if mode == 'u':
val = (amplitude * nm.sin(time * omega) * nm.sin(k * y + shift))

elif mode == 'du':
val = (amplitude * omega * nm.cos(time * omega)
* nm.sin(k * y + shift))

elif mode == 'ddu':
val = (-amplitude * omega**2 * nm.sin(time * omega)
* nm.sin(k * y + shift))

return val

functions = {
'get_u' : (lambda ts, coor, **kwargs: get_ebcs(ts, coor),),
'get_du' : (lambda ts, coor, **kwargs: get_ebcs(ts, coor, mode='du'),),
'get_ddu' : (lambda ts, coor, **kwargs: get_ebcs(ts, coor, mode='ddu'),),
}

ebcs = {
'Seismic' : ('Bottom', {'u.2' : 'get_u', 'du.2' : 'get_du',
'ddu.2' : 'get_ddu'}),
'Pot0' : ('Bottom', {'p.all' : 0.0}),
}

ics = {
'ic' : ('Omega', {'u.all' : 0.0, 'du.all' : 0.0, 'p.0' : 0.0}),
}

equations = {
'1' : """dw_dot.i.Omega(m.rho, v, ddu)
+ dw_lin_elastic.i.Omega(m.C, v, u)
- dw_piezo_coupling.i.Omega(m.e, v, p)
= 0""",
'2' : """dw_piezo_coupling.i.Omega(m.e, u, q)
+ dw_diffusion.i.Omega(m.kappa, q, p)
- de_surface_flux.i.Top(m.kappa, q, p)
+ de_surface_flux.i.Top(m.kappa, p, q)
- de_surface_flux.i.Top(m.kappa, pc, q)
+ dw_dot.i.Top(m.penalty, q, p)
- dw_dot.i.Top(m.penalty, q, pc)
= 0""",
'3' : """de_surface_flux.i.Top(m.kappa, qc, dp/dt)
+ 0.5 * dw_dot.i.Top(m.iR, qc, pc)
+ 0.5 * dw_dot.i.Top(m.iR, qc, pc[-1])
= 0""",
}

solvers = {
'lsd' : ('ls.auto_direct', {
# Reuse the factorized linear system from the first time step.
'use_presolve' : True,
# Speed up the above by omitting the matrix digest check used
# normally for verification that the current matrix corresponds to
# the factorized matrix stored in the solver instance. Use with
# care!
'use_mtx_digest' : False,
# Increase when getting MUMPS error -9.
'memory_relaxation' : 20,
}),
'newton' : ('nls.newton', {
'i_max'      : 1,
'eps_a'      : 1e-6,
'eps_r'      : 1e-6,
'ls_on'      : 1e100,
}),
'tsn' : ('ts.newmark', {
't0' : 0.0,
't1' : t1,
'dt' : dt,
'n_step' : None,

'is_linear' : True,
# Without this the adaptive time-stepping cannot work.

'beta' : 0.25,
'gamma' : 0.5,

'verbose' : 1,
}),
'tscedl' : ('tsc.ed_linear', {
'eps_r' : (1e-4, 1e-2),
'eps_a' : (1e-8, 1e-3),
'fmin' : 0.3,
'fmax' : 2.5,
'fsafety' : 0.85,
'red_factor' : 0.9,
'inc_wait' : 10,
'min_inc_factor' : 1.5,
}),
}

options = {
'ts' : tss_name,
'tsc' : tsc_name if adaptive else None,
'nls' : 'newton',
'ls' : ls_name,

'save_times' : save_times,

'active_only' : active_only,
'auto_transform_equations' : True,

'output_format' : 'h5',
'output_dir' : output_dir,
'post_process_hook' : _post_process,
'post_process_hook_final' : _plot_voltage,
}

return locals()