# Electromechanical Applications

(Difference between revisions)
 Revision as of 14:08, 28 April 2011 (view source)Amir (Talk | contribs) (→Theory)← Older edit Latest revision as of 20:02, 28 April 2011 (view source)Amir (Talk | contribs) (→Weak form) (42 intermediate revisions by one user not shown) Line 1: Line 1: == General Description == == General Description == − We have provided the electro-mechanical applications in [[Kratos | KRATOS]] for the simulations of [http://en.wikipedia.org/wiki/Ferroelectricity  ferroelectric] and [http://en.wikipedia.org/wiki/Piezoelectricity  piezoelectric] materials. The code has been developed mainly to study the crack propagation in these materials but it can be effectively used for other purposes. Here we start from a brief theory of these materials and then we move to the numerical implementation in [[Kratos | KRATOS]]. The theory of ferroelectric materials is presented in our published paper which can be accessed in the following link. + We have provided the electro-mechanical applications in [[Kratos | KRATOS]] for the simulations of [http://en.wikipedia.org/wiki/Ferroelectricity  ferroelectric] and [http://en.wikipedia.org/wiki/Piezoelectricity  piezoelectric] materials. The code has been developed mainly to study the crack propagation in these materials but it can be effectively used for other purposes regarding an intact model. An excellent forum regarding the fracture of ferroelectrics is presented in [http://imechanica.org/node/3064 imechanica - Fracture of Ferroelectrics]. Here we start from a brief theory of these materials and then we move to the numerical implementation in [[Kratos | KRATOS]]. The theory of ferroelectric materials is presented in our published paper which can be accessed through the following link. http://dx.doi.org/10.1016/j.actamat.2011.03.030 http://dx.doi.org/10.1016/j.actamat.2011.03.030 Line 8: Line 8: ferroelectric single crystals. Acta Mater (2011), doi:10.1016/j.actamat.2011.03.030 ferroelectric single crystals. Acta Mater (2011), doi:10.1016/j.actamat.2011.03.030 − === Theory === − The total electro-mechanical enthalpy for a ferroelectric body occupying a region is stated as: + All the theory and the developed codes for the numerical implementation are based on the model presented in this paper. + + + === Theory === + ==== Ferroelectric Materials ==== + The total electro-mechanical enthalpy for a ferroelectric body occupying a region $\Omega$ is stated as: $[itex] Line 34: Line 38: [itex] [itex] − U(p_{i,j})= \frac{a_0}{2}(p^2_{1,1} + p^2_{1,2} + p^2_{2,1} + p^2_{2,2}) + U(p_{i,j})= \frac{a_0}{2}(p^2_{1,1} + p^2_{1,2} + p^2_{2,1} + p^2_{2,2}), \quad \quad \quad \quad \quad (4)$ [/itex] Line 42: Line 46: W(p_i, \varepsilon_{jk})= -\frac{b_1}{2}(\varepsilon_{11}p^2_{1} + \varepsilon_{22}p^2_{2}) - \frac{b_2}{2}(\varepsilon_{11}p^2_{2} + \varepsilon_{22}p^2_{1}) - {b_3}(\varepsilon_{21} + \varepsilon_{12})p_{1}p_{2} W(p_i, \varepsilon_{jk})= -\frac{b_1}{2}(\varepsilon_{11}p^2_{1} + \varepsilon_{22}p^2_{2}) - \frac{b_2}{2}(\varepsilon_{11}p^2_{2} + \varepsilon_{22}p^2_{1}) - {b_3}(\varepsilon_{21} + \varepsilon_{12})p_{1}p_{2} +\frac{c_1}{2}(\varepsilon^2_{11} + \varepsilon^2_{22}) + {c_2}\varepsilon_{11}\varepsilon_{22} + \frac{c_3}{2}(\varepsilon^2_{12} + \varepsilon^2_{21}), +\frac{c_1}{2}(\varepsilon^2_{11} + \varepsilon^2_{22}) + {c_2}\varepsilon_{11}\varepsilon_{22} + \frac{c_3}{2}(\varepsilon^2_{12} + \varepsilon^2_{21}), − [/itex] + \quad  \quad  \quad  \quad \quad (5) + [/itex] Line 48: Line 53: $[itex] \chi(p_i) = \alpha_1(p^2_{1} + p^2_{2}) + \alpha_{11}(p^4_{1} + p^4_{2}) + \alpha_{12}(p^2_{1}p^2_{2}) + \alpha_{111}(p^6_{1} + p^6_{2}) + \alpha_{112}(p^2_{1}p^4_{2} + p^2_{2}p^4_{1}) \chi(p_i) = \alpha_1(p^2_{1} + p^2_{2}) + \alpha_{11}(p^4_{1} + p^4_{2}) + \alpha_{12}(p^2_{1}p^2_{2}) + \alpha_{111}(p^6_{1} + p^6_{2}) + \alpha_{112}(p^2_{1}p^4_{2} + p^2_{2}p^4_{1}) − + \alpha_{1111}(p^8_{1} + p^8_{2}) + \alpha_{1112}(p^6_{1}p^2_{2} + p^6_{2}p^2_{1}) + \alpha_{1122}(p^4_{1}p^4_{2}), + + \alpha_{1111}(p^8_{1} + p^8_{2}) + \alpha_{1112}(p^6_{1}p^2_{2} + p^6_{2}p^2_{1}) + \alpha_{1122}(p^4_{1}p^4_{2}), \quad (6)$ [/itex] − where the combination of energy functions $\chi$ and $W$ is the + The combination of energy functions $\chi$ and $W$ is the total Landau-Devonshire energy density, $a_0$ is the scaling total Landau-Devonshire energy density, $a_0$ is the scaling parameter of the domain wall energy, $b_i (i=1,2,3)$ are the parameter of the domain wall energy, $b_i (i=1,2,3)$ are the − constants of the coupling terms between strain and polarization and + constants of the coupling terms between strain and polarization, − $c_i (i=1,2,3)$ are the elastic constants. + $c_i (i=1,2,3)$ are the elastic constants and $\alpha_{1}$, $\alpha_{11}$, $\alpha_{12}$, $\alpha_{111}$, $\alpha_{112}$, $\alpha_{1111}$, $\alpha_{1112}$, $\alpha_{1122}$, are the constants of the phase-separation potential $\chi$. + + This model leads to six degrees of freedom  $\mathbf{u}$, $\mathbf{p}$, $\phi$ and $v$ per node in the case of plane polarization and strain. + + ==== Weak form ==== + + In this model $v$ is selected together with the polarization $\mathbf{p}$ as primary order parameters and finite mobilities are introduced for the micro-structure and fracture processes. Thus, the weak form of the gradient flow for the primary variables, together with the equations for mechanical and electrostatic equilibria follow from + + $+ \mu_p \int_\Omega \dot{p_i}\delta p_i ~\mathrm{d}\Omega = - \delta H[\mathbf{u},v,\mathbf{p},\phi; \delta \mathbf{p}] + = - \int_\Omega \frac{\partial h}{\partial p_i}\delta p_i ~\mathrm{d}\Omega, \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad (7) +$ + + $+ \mu_v \int_\Omega \dot{v}\delta v ~\mathrm{d}\Omega = - \delta H[\mathbf{u},v,\mathbf{p},\phi; \delta v] = - \int_\Omega \frac{\partial h}{\partial v}\delta v ~\mathrm{d}\Omega - 2G_c \int_\Omega \left(\frac{v-1}{4\kappa} \delta v + \kappa{v_{,i}}\delta{v_{,i}}\right) ~\mathrm{d}\Omega, \quad \quad \quad \quad\quad \quad \quad (8) +$ + + $+ 0 = \delta H[\mathbf{u},v,\mathbf{p},\phi; \delta \mathbf{u}] = \int_\Omega \frac{\partial h}{\partial \varepsilon_{ij}}\delta\varepsilon_{ij} ~\mathrm{d}\Omega - \int_{\Gamma_{N,\mathbf{u}}} t_i\delta u_i ~\mathrm{d}S, \quad \quad ----> \text{Mechanical Equilibrium} \quad\quad\quad\quad (9) +$ + + $+ 0 = -\delta H[\mathbf{u},v,\mathbf{p},\phi; \delta \phi] = -\int_\Omega \frac{\partial h}{\partial E_i} \delta E_{i} ~\mathrm{d}\Omega - \int_{\Gamma_{N,\phi}} \omega\delta \phi ~\mathrm{d}S, \quad ----> \text{Electrostatic Equilibrium} \quad\quad\quad (10) +$ + + where $1/{\mu_p} > 0$ and $1/{\mu_v} > 0$ are the mobilities of the + processes. The weak form of the evolution and equilibrium + equations is discretized in space with standard finite elements. + Equations (7) and (8) are discretized in time + with a semi-implicit scheme. + + '''NOTE''' that  for the simulations of ferroelectric materials without considering a crack, the scalar field $v$ should be fixed to 1 in all the domain. In this case, the degrees of freedom reduces to five and there is no need to solve equation (8). + + A simple algorithm to solve the coupled system in a straightforward staggered approach is presented in [http://dx.doi.org/10.1016/j.actamat.2011.03.030]. This algorithm is based on solving equations (7) - (10) in an iterative process which leads to the coupled microstructure and fracture evolution. We implement this algorithm in a [[Python]] interface file which is presented in the next section. + + === Numerical Implementation === + + In this section, we present the structure of the electromechanical element in [[Kratos | KRATOS]] where each part of the code is discussed according to  the theory in the previous section. The Python run file is also presented in this section. − This model leads to to six degrees of freedom  $\mathbf{u}$, $\mathbf{p}$, $\phi$ and $v$ per node in the case of plane polarization and strain. NOTE that for the simulations of ferroelectric materials without considering a crack, the scalar field $v$ should be fixed to 1 in all the domain. In this case, the degrees of freedom reduces to five. + ==== Electromechanical Element: ElectroMechanical2D.cpp ====

## General Description

We have provided the electro-mechanical applications in KRATOS for the simulations of ferroelectric and piezoelectric materials. The code has been developed mainly to study the crack propagation in these materials but it can be effectively used for other purposes regarding an intact model. An excellent forum regarding the fracture of ferroelectrics is presented in imechanica - Fracture of Ferroelectrics. Here we start from a brief theory of these materials and then we move to the numerical implementation in KRATOS. The theory of ferroelectric materials is presented in our published paper which can be accessed through the following link.

Abdollahi A, Arias I. Phase-field modeling of the coupled microstructure and fracture evolution in ferroelectric single crystals. Acta Mater (2011), doi:10.1016/j.actamat.2011.03.030

All the theory and the developed codes for the numerical implementation are based on the model presented in this paper.

### Theory

#### Ferroelectric Materials

The total electro-mechanical enthalpy for a ferroelectric body occupying a region Ω is stated as:

$H[\mathbf{u},v,\mathbf{p},\phi] = \int_\Omega h(\mathbf{\varepsilon}(\mathbf{u}),\mathbf{p},\nabla\mathbf{p},\mathbf{E}(\phi),v)~\mathrm{d}\Omega + G_c\int_{\Omega} \left[\frac{(1-v)^2}{4\kappa} + \kappa|\nabla{v}|^2\right]~\mathrm{d}\Omega - \int_{\Gamma_{N,\mathbf{u}}} \mathbf{t}\cdot\mathbf{u}~{\rm d}S + \int_{\Gamma_{N,\phi}} \omega\phi~{\rm d}S \quad \quad \quad \quad \quad (1)$

where $\mathbf{t}$ and ω are the tractions and surface charge density respectively, and $\Gamma_{N,\mathbf{u}}$ and ΓN are the parts of the boundary of the domain $\partial\Omega$ where mechanical and electrical Neumann boundary conditions are applied. $\mathbf{\varepsilon}$ is the strain tensor associated with the mechanical displacement $\mathbf{u}$, $\mathbf{\varepsilon} = 1/2(\nabla\mathbf{u} + \nabla^{T}\mathbf{u})$, $\mathbf{p}$ is the polarization, $\mathbf{E}$ is the electric field defined as $\mathbf{E} = -\nabla\phi$, where φ is the electric potential. Gc is the critical energy release rate or the surface energy density in Griffith's theory. The scalar field v provides a diffuse representation of the fracture zone, κ is a positive regularization constant to regulate the size of the fracture zone. The electro-mechanical enthalpy density h considering permeable and impermeable cracks follows

$h(\mathbf{\varepsilon},\mathbf{p},\nabla\mathbf{p},\mathbf{E},v) = (v^2 + \eta_\kappa)\left[U(\nabla \mathbf{p}) + W(\mathbf{p},\mathbf{\varepsilon})\right] + \chi(\mathbf{p}) - \frac{\varepsilon_0}{2}|\mathbf{E}|^2 -\mathbf{E}\cdot\mathbf{p} \quad \quad \quad \quad \quad \quad \quad \quad (2) \quad \text{Permeable Crack}$

$h(\mathbf{\varepsilon},\mathbf{p},\nabla\mathbf{p},\mathbf{E},v) = (v^2 + \eta_\kappa)\left[U(\nabla \mathbf{p}) + W(\mathbf{p},\mathbf{\varepsilon}) - \frac{\varepsilon_0}{2}|\mathbf{E}|^2 -\mathbf{E}\cdot\mathbf{p}\right] + \chi(\mathbf{p}) \quad \quad \quad \quad \quad \quad (3) \quad \text{Impermeable Crack}$

where the energy functions U, W and χ are stated as

$U(p_{i,j})= \frac{a_0}{2}(p^2_{1,1} + p^2_{1,2} + p^2_{2,1} + p^2_{2,2}), \quad \quad \quad \quad \quad (4)$

$W(p_i, \varepsilon_{jk})= -\frac{b_1}{2}(\varepsilon_{11}p^2_{1} + \varepsilon_{22}p^2_{2}) - \frac{b_2}{2}(\varepsilon_{11}p^2_{2} + \varepsilon_{22}p^2_{1}) - {b_3}(\varepsilon_{21} + \varepsilon_{12})p_{1}p_{2} +\frac{c_1}{2}(\varepsilon^2_{11} + \varepsilon^2_{22}) + {c_2}\varepsilon_{11}\varepsilon_{22} + \frac{c_3}{2}(\varepsilon^2_{12} + \varepsilon^2_{21}), \quad \quad \quad \quad \quad (5)$

$\chi(p_i) = \alpha_1(p^2_{1} + p^2_{2}) + \alpha_{11}(p^4_{1} + p^4_{2}) + \alpha_{12}(p^2_{1}p^2_{2}) + \alpha_{111}(p^6_{1} + p^6_{2}) + \alpha_{112}(p^2_{1}p^4_{2} + p^2_{2}p^4_{1}) + \alpha_{1111}(p^8_{1} + p^8_{2}) + \alpha_{1112}(p^6_{1}p^2_{2} + p^6_{2}p^2_{1}) + \alpha_{1122}(p^4_{1}p^4_{2}), \quad (6)$

The combination of energy functions χ and W is the total Landau-Devonshire energy density, a0 is the scaling parameter of the domain wall energy, bi(i = 1,2,3) are the constants of the coupling terms between strain and polarization, ci(i = 1,2,3) are the elastic constants and α1, α11, α12, α111, α112, α1111, α1112, α1122, are the constants of the phase-separation potential χ.

This model leads to six degrees of freedom $\mathbf{u}$, $\mathbf{p}$, φ and v per node in the case of plane polarization and strain.

#### Weak form

In this model v is selected together with the polarization $\mathbf{p}$ as primary order parameters and finite mobilities are introduced for the micro-structure and fracture processes. Thus, the weak form of the gradient flow for the primary variables, together with the equations for mechanical and electrostatic equilibria follow from

$\mu_p \int_\Omega \dot{p_i}\delta p_i ~\mathrm{d}\Omega = - \delta H[\mathbf{u},v,\mathbf{p},\phi; \delta \mathbf{p}] = - \int_\Omega \frac{\partial h}{\partial p_i}\delta p_i ~\mathrm{d}\Omega, \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad (7)$

$\mu_v \int_\Omega \dot{v}\delta v ~\mathrm{d}\Omega = - \delta H[\mathbf{u},v,\mathbf{p},\phi; \delta v] = - \int_\Omega \frac{\partial h}{\partial v}\delta v ~\mathrm{d}\Omega - 2G_c \int_\Omega \left(\frac{v-1}{4\kappa} \delta v + \kappa{v_{,i}}\delta{v_{,i}}\right) ~\mathrm{d}\Omega, \quad \quad \quad \quad\quad \quad \quad (8)$

$0 = \delta H[\mathbf{u},v,\mathbf{p},\phi; \delta \mathbf{u}] = \int_\Omega \frac{\partial h}{\partial \varepsilon_{ij}}\delta\varepsilon_{ij} ~\mathrm{d}\Omega - \int_{\Gamma_{N,\mathbf{u}}} t_i\delta u_i ~\mathrm{d}S, \quad \quad ----> \text{Mechanical Equilibrium} \quad\quad\quad\quad (9)$

$0 = -\delta H[\mathbf{u},v,\mathbf{p},\phi; \delta \phi] = -\int_\Omega \frac{\partial h}{\partial E_i} \delta E_{i} ~\mathrm{d}\Omega - \int_{\Gamma_{N,\phi}} \omega\delta \phi ~\mathrm{d}S, \quad ----> \text{Electrostatic Equilibrium} \quad\quad\quad (10)$

where 1 / μp > 0 and 1 / μv > 0 are the mobilities of the processes. The weak form of the evolution and equilibrium equations is discretized in space with standard finite elements. Equations (7) and (8) are discretized in time with a semi-implicit scheme.

NOTE that for the simulations of ferroelectric materials without considering a crack, the scalar field v should be fixed to 1 in all the domain. In this case, the degrees of freedom reduces to five and there is no need to solve equation (8).

A simple algorithm to solve the coupled system in a straightforward staggered approach is presented in [1]. This algorithm is based on solving equations (7) - (10) in an iterative process which leads to the coupled microstructure and fracture evolution. We implement this algorithm in a Python interface file which is presented in the next section.

### Numerical Implementation

In this section, we present the structure of the electromechanical element in KRATOS where each part of the code is discussed according to the theory in the previous section. The Python run file is also presented in this section.