Voigt notations and electrical potential boundaries for linear piezoelectricity

[Ycliu6]

Hello @kinnala,

I hope you are doing well!

I am a novice in either fem or skfem, and currently working on solving eigenvalues for linear piezoelectricity. The governing equations are presented:
$\nabla\cdot(C\varepsilon(u)+e^T\nabla(\phi)) = \rho\omega^2u$ in $\Omega$
$\nabla\cdot(e\varepsilon(u)-\kappa\nabla(\phi)) = 0$ in $\Omega$
where $\varepsilon(*)$ is sym_grad $(*)$ denotes strain, $\phi$ is electric potential.
If ignore the traction on the surface first, their weak form is like:
$\int_\Omega(C\varepsilon(u)): \varepsilon(p) dx$ + $\int_\Omega(e^T\nabla\phi): \varepsilon(p) dx$ = $\int_\Omega\rho\omega^2 u \cdot p dx$
$\int_\Omega(e\varepsilon(u)): \nabla q dx$ - $\int_\Omega(\kappa \nabla \phi): \nabla q dx$ = 0
where $u$, $\phi$ are trial functions and $p$, $q$ are test functions.
For further details, please refer to: https://doi.org/10.1093/oso/9780198864721.003.0018 (p336)

Could I kindly ask a few questions,
Q1: Could you please suggest the proper application of the Voigt notation to BilinearForm. I have attempted to use the code below based on a beam model:
`from skfem import *
import numpy as np
import matplotlib.pyplot as plt
import skfem.visuals.matplotlib
from skfem.helpers import ddot, sym_grad, grad, dot, mul

mesh = MeshHex.init_tensor(
np.linspace(0, 1, 5),
np.linspace(-.05, .05, 5),
np.linspace(-.05, .05, 5),
).with_boundaries({
'left': lambda x: x[0] == 0.,
'right': lambda x: x[0] == 1.,
})
uelem = ElementVector(ElementHex1()) # Element for displacement
phielem = ElementHex1() # Element for electric potential
elems = {"u": uelem, "phi": phielem}
basis = {
field: Basis(mesh, e, intorder=2)
for field, e in elems.items()
}

#Define material properties
stiffness = np.array([[1.68001563e+11, 1.10351384e+11, 9.99053252e+10, 0, 0, 0],
[1.10351384e+11, 1.68001563e+11, 9.99053252e+10, 0, 0, 0],
[9.99053252e+10, 9.99053252e+10, 1.22636143e+11, 0, 0, 0],
[0, 0, 0, 3.01260000e+10, 0, 0],
[0, 0, 0, 0, 3.01260000e+10, 0],
[0, 0, 0, 0, 0, 2.88251724e+10]])
piezoelectric = np.array([[0, 0, 0, 0, 9.8568, 0],
[0, 0, 0, 9.8568, 0, 0],
[-2.80119191, -2.80119191, 14.6914, 0, 0, 0]])
permittivity = np.array([[1.05634002e-08, 0, 0],
[0, 1.05634002e-08, 0],
[0, 0, 1.14396107e-08]])
density = 7700

#tensor (3x3) to vector (6x1)
def ten2vec(T):
return np.array([T[0,0], T[1,1], T[2,2], 2T[1,2], 2T[0,2], 2*T[0,1]])
#vector (6x1) to tenstor (3x3)
def vec2ten(T):
return np.array([[T[0], T[5], T[4]],[T[5], T[1], T[3]],[T[4], T[3], T[2]]])

@BilinearForm
def K00(u, v, w):
vec_symu = ten2vec(sym_grad(u))
stif_mul = mul(stiffness, vec_symu)
ten_symu = vec2ten(stif_mul)
return ddot(ten_symu, sym_grad(v))

@BilinearForm
def K10(u, v, w):
vec_symu = ten2vec(sym_grad(u))
piezo_mul = mul(piezoelectric, vec_symu)
return ddot(piezo_mul, grad(v))

@BilinearForm
def K11(u, v, w):
return ddot(mul(permittivity, grad(u)), grad(v))

#Stiffness matrix
K00 = asm(K00, basis['u'], basis['u'])
K10 = asm(K10, basis['u'], basis['phi'])
K11 = asm(K11, basis['phi'],basis['phi'])
K = bmat(
[[K00, K10.T],
[K10, -K11]], "csr"
)

@BilinearForm
def mass_u(u, v, w):
return dot(density * u, v)

@BilinearForm
def mass_phi(u, v, w):
return 0 * u * v

#Mass matrix
M_u = asm(mass_u, basis['u'], basis['u'])
M_phi = asm(mass_phi, basis['phi'], basis['phi'])
M = bmat(
[[M_u, None],
[None, M_phi]], 'csr'
)
#Free boundaries for u
L, x = solve(K, M, solver=solver_eigen_scipy_sym(k=50, sigma=0.0))
`

Q2: I am wondering how to apply electrical potential on the left and right facets for the problem above to form an electric field. Should I be using LinearForm as ex18 or ex36? If so, do I need to create a facet basis?

Thanks so much for your help in advance, and look forward to hearing from you.

Cheers,

Yuchen Liu

[kinnala]

Q1. What you describe should work. However, I'd probably implement something like voigtmul which does this tensor-2-vector mapping, matrix multiplication and reverse mapping in one step, and then pretend all those tensors are written in Voigt notation. I personally do not use Voigt notation that much, but I guess if your reference is written that way then you are inclined to follow it also.

[Ycliu6]

Thanks so much for your response, really appreciate it. That would be a nice way to implement 'voigtmul' as you mentioned! And no rush for the Q2.