diff --git a/2d_polarization.py b/2d_polarization.py new file mode 100644 index 0000000000000000000000000000000000000000..18c1fd7519b12996f7dc45072a86c9e80fa6f4fc --- /dev/null +++ b/2d_polarization.py @@ -0,0 +1,207 @@ +""" +Material: BTO. +Order Parameter: Spontaneous Polarization. +Model: Phase Field with external electric field and strain. +Numerics: Fourier spectral method. 2D. + +You can store all data as .h5 (Write_to_HDF5) + +""" + +from material_data import Mat_Data_BTO_D, initial_polarization +from fourier_frequencies import fourier_space +from solve_piezoelectricity import Solve_Piezoelectricity_2D +from green_operator import Green_Operator_Piezoelectricity_2D +from calculate_derivatives import Calculate_Derivatives_2D, Calculate_Derivatives_Land_Potential_2D +from make_plot import Plot_Polarization, Plot_Potential +from store_hdf5 import Write_to_HDF5 +# ---------------------------------------------------------------------------- +import numpy as np +import timeit +import itertools +#----------------------------------------------------------------------------- + +def Evolve_Sponteneous_Polarization_2D(folder, Nx, Ny, nsteps, nt, dt, domain_type, E_app_y, E_app_x, eps_app_x, eps_app_y): + + # ---------------------------------------------------------------------------- + # import material data + G, l_wall, P0, mob, e00, k_sep, k_grad, l_coef, e_piezo, C0, C_cub, K0 = Mat_Data_BTO_D() + dx, dy = (l_wall/5, l_wall/5) + # ---------------------------------------------------------------------------- + # initial polarization + P = initial_polarization(Nx, Ny, domain_type) * P0 + # ---------------------------------------------------------------------------- + # frequencies + freq = fourier_space(Nx, dx, Ny, dy) + # ---------------------------------------------------------------------------- + # save inital polarization + Plot_Polarization(folder + "/P_", 0, Nx, Ny, P[0], P[1], P0) + #Write_to_HDF5(folder, 0, P) + # ---------------------------------------------------------------------------- + # plot the potential + Plot_Potential(folder, P0, Nx, l_coef[0], l_coef[1], l_coef[2], l_coef[3], l_coef[4]) + # ---------------------------------------------------------------------------- + # assign tensor for each grid + C0 = np.einsum('ijkl, xy -> ijklxy', C0, np.ones((Nx,Ny))) + K0 = np.einsum('ij, xy -> ijxy', K0, np.ones((Nx,Ny))) + # ---------------------------------------------------------------------------- + eps_ext = np.zeros((2,2,Nx,Ny)) # applied elas. field + E_ext = np.zeros((2,Nx,Ny)) # applied elec. field + # ---------------------------------------------------------------------------- + eps_ext[0,0] = eps_app_x; eps_ext[1,1] = eps_app_y + E_ext[0] = E_app_x; E_ext[1] = E_app_x + # ---------------------------------------------------------------------------- + # random electric field, zero for perfect crystal + E_random = np.zeros_like(E_ext) + # ---------------------------------------------------------------------------- + # derivatives + eP_val, deP_dPx_val, deP_dPy_val, eps0_val, deps0_P_dPx_val, deps0_P_dPy_val = Calculate_Derivatives_2D() + Psi_val, dPsi_dPx_val, dPsi_dPy_val = Calculate_Derivatives_Land_Potential_2D() + # ---------------------------------------------------------------------------- + # time evolution + for step in range(nsteps): + # ------------------------------------------------------------------------ + print('\n---------------------- Time step:\t' + str(step)+'\t----------------------') + # ------------------------------------------------------------------------- + # calculate derivatives, convert symbolic derivatives into numpy + eP = np.array(eP_val(P[0], P[1], P0, e_piezo[0], e_piezo[1], e_piezo[2])) + deP_dPx = np.array(deP_dPx_val(P[0], P[1], P0, e_piezo[0], e_piezo[1], e_piezo[2])) + deP_dPy = np.array(deP_dPy_val(P[0], P[1], P0, e_piezo[0], e_piezo[1], e_piezo[2])) + # ------------------------------------------------------------------------- + eps0 = np.array(eps0_val(P[0], P[1], P0, e00)) + deps0_P_dPx = np.array(deps0_P_dPx_val(P[0], P[1], P0, e00)) + deps0_P_dPy = np.array(deps0_P_dPy_val(P[0], P[1], P0, e00)) + # ------------------------------------------------------------------------- + e0 = np.zeros((2,2,2,Nx,Ny)) # homogeneous piez. tensor, volume average of e(P) + for i,j,k in itertools.product(range(2), repeat=3): + e0[i,j,k] = np.mean(eP[i,j,k]) + # ------------------------------------------------------------------------- + # greens functions + G_elas, G_piezo, G_dielec = Green_Operator_Piezoelectricity_2D(C0, e0, K0, freq) + # ------------------------------------------------------------------------ + # Solve piezoelectricity + sigma, D, eps_elas, E = Solve_Piezoelectricity_2D(C0, K0, e0, C0, K0, eP, G_elas, G_piezo, G_dielec, + eps0, eps_ext, E_ext, E_random, Nx, Ny, P) + # ------------------------------------------------------------------------ + del G_elas, G_piezo, G_dielec + # Print energies + # ------------------------------------------------------------------------ + H_elas = 0.5*np.einsum('ijkl...,ij...,kl...',C0, eps_elas, eps_elas) + # ------------------------------------------------------------------------ + print('\nElastic Energy:\t\t\t\t\t\t'+ str(np.mean(H_elas)) + '\t[N/m2]') + # ------------------------------------------------------------------------ + H_piezo = -1*np.einsum('ijk...,ij...,k...', eP, eps_elas, E) + # ------------------------------------------------------------------------ + print('Piezoelectric Energy:\t\t\t\t\t'+ str(np.mean(H_piezo)) + '\t[N/m2]') + # ------------------------------------------------------------------------ + H_elec = - 0.5*np.einsum('ij...,i...,j...', K0, E, E) - np.einsum('i...,i...', P, E) + # ------------------------------------------------------------------------ + print('Electric Energy:\t\t\t\t\t'+ str(np.mean(H_elec)) + '\t[N/m2]') + # ------------------------------------------------------------------------ + # Derivatives + deps_elas_dPx = -1*deps0_P_dPx; deps_elas_dPy = -1*deps0_P_dPy + # ------------------------------------------------------------------------ + dH_elas_dPx = 0.5*(np.einsum('ijkl...,ij...,kl...', C0, deps_elas_dPx, eps_elas) + + np.einsum('ijkl...,ij...,kl...', C0, eps_elas, deps_elas_dPx) ) + dH_elas_dPy = 0.5*(np.einsum('ijkl...,ij...,kl...', C0, deps_elas_dPy, eps_elas) + + np.einsum('ijkl...,ij...,kl...', C0, eps_elas, deps_elas_dPy) ) + # ------------------------------------------------------------------------ + dH_piezo_dPx = -1*(np.einsum('ijk...,ij...,k...', deP_dPx, eps_elas, E) + + np.einsum('ijk...,ij...,k...', eP, deps_elas_dPx, E) ) + dH_piezo_dPy = -1*(np.einsum('ijk...,ij...,k...', deP_dPy, eps_elas, E) + + np.einsum('ijk...,ij...,k...', eP, deps_elas_dPy, E) ) + # ------------------------------------------------------------------------ + dH_elec_dPx = -1*E[0] + dH_elec_dPy = -1*E[1] + # ------------------------------------------------------------------------ + Psi = np.array(Psi_val(P[0], P[1], P0, l_coef[0], l_coef[1], l_coef[2], l_coef[3], l_coef[4])) + # ------------------------------------------------------------------------ + H_sep = k_sep * G/l_wall * Psi + # ------------------------------------------------------------------------ + print('Separation Energy:\t\t\t\t\t'+ str(np.mean(H_sep)) + '\t[N/m2]') + # ------------------------------------------------------------------------ + # ------------------------------------------------------------------------ + dPsi_dPx = np.array(dPsi_dPx_val(P[0], P[1], P0, l_coef[1], l_coef[2], l_coef[3], l_coef[4])) + dPsi_dPy = np.array(dPsi_dPy_val(P[0], P[1], P0, l_coef[1], l_coef[2], l_coef[3], l_coef[4])) + # ------------------------------------------------------------------------ + dH_dPx = k_sep * G/l_wall * dPsi_dPx + dH_elas_dPx + dH_piezo_dPx + dH_elec_dPx + dH_dPy = k_sep * G/l_wall * dPsi_dPy + dH_elas_dPy + dH_piezo_dPy + dH_elec_dPy + # ------------------------------------------------------------------------ + dH_dPx_k = np.fft.fft2(dH_dPx) + dH_dPy_k = np.fft.fft2(dH_dPy) + # ------------------------------------------------------------------------ + Px_k = np.fft.fft2(P[0]); Py_k = np.fft.fft2(P[1]) + # ------------------------------------------------------------------------ + dPx_dx = np.fft.ifft2(1j*freq[0]*Px_k).real + dPx_dy = np.fft.ifft2(1j*freq[1]*Px_k).real + dPy_dx = np.fft.ifft2(1j*freq[0]*Py_k).real + dPy_dy = np.fft.ifft2(1j*freq[1]*Py_k).real + # ------------------------------------------------------------------------ + H_grad = 0.5*k_grad*G*l_wall/P0**2*(dPx_dx**2 + dPx_dy**2 + dPy_dx**2 + dPy_dy**2) + print('Gradient Energy:\t\t\t\t'+str(np.mean(H_grad))) + # ------------------------------------------------------------------------ + del dH_elas_dPx, dH_elas_dPy, dH_piezo_dPx, dH_piezo_dPy, dH_elec_dPx, dH_elec_dPy + del dH_dPx, dH_dPy, dPx_dx, dPx_dy, dPy_dx, dPy_dy, dPsi_dPx, dPsi_dPy + # ------------------------------------------------------------------------ + Px_k = (Px_k - dt * mob * dH_dPx_k) / (1 + dt * mob * k_grad *G*l_wall/P0**2 * (freq[0]**2 + freq[1]**2)) + Py_k = (Py_k - dt * mob * dH_dPy_k) / (1 + dt * mob * k_grad *G*l_wall/P0**2 * (freq[0]**2 + freq[1]**2)) + # ------------------------------------------------------------------------ + Px = np.fft.ifft2(Px_k).real + Py = np.fft.ifft2(Py_k).real + # ------------------------------------------------------------------------ + P = np.array([Px, Py]) + # ------------------------------------------------------------------------ + if ((step + 1) % nt == 0 ): + # -------------------------------------------------------------------- + Plot_Polarization(folder + "/P_", step + 1, Nx, Ny, P[0], P[1], P0) + """ + Use the following if you want to save any of them: + Stress=sigma + Elas_strain=eps_elas + Elec_disp=D + Elec_field=E + Elas_en=H_elas + Grad_en=H_grad + Sep_en=H_sep + Elec_en=H_elec + Piezo_en=H_piezo + """ + #Write_to_HDF5(folder, step + 1, P)#, Grad_en=H_grad) + # -------------------------------------------------------------------- + del sigma, D, eps_elas, E + del H_elas, H_piezo, H_elec, H_grad, Psi, H_sep + # ---------------------------------------------------------------------------- + +def main(): + + start_tm = timeit.default_timer() + # ---------------------------------------------------------------------------- + # grid parameters + Nx, Ny = (360, 360) + # time step parameters + nsteps, nt, dt = (50000, 100, 1e-12) # dt is in sec + # applied fields in each direction + E_app_x, E_app_y = (0, 0) # must be in V/m + eps_app_x, eps_app_y = (0, 0) + # ---------------------------------------------------------------------------- + # folder to save results + import os + folder = os.getcwd() + # ---------------------------------------------------------------------------- + domain_type = "random" # domain_type = "ramdom" OR "90" OR "180" + # ---------------------------------------------------------------------------- + # Evolve + Evolve_Sponteneous_Polarization_2D(folder, Nx, Ny, nsteps, nt, dt, domain_type, E_app_y, E_app_x, eps_app_x, eps_app_y) + # ---------------------------------------------------------------------------- + stop_tm = timeit.default_timer() + print('\nExecution time: ' + str( format((stop_tm - start_tm)/60,'.5f')) + ' min\n' ) + + +if __name__ == "__main__": + main() + + + + + + \ No newline at end of file diff --git a/calculate_derivatives.py b/calculate_derivatives.py new file mode 100644 index 0000000000000000000000000000000000000000..a523485b337cce26a854b1587d567960b381f032 --- /dev/null +++ b/calculate_derivatives.py @@ -0,0 +1,93 @@ +""" +Symbolic differentiation of the spontaneous strain and piezoelectric tensor wrt P +""" + +import sympy as sym +import numpy as np +import itertools + +def Calculate_Derivatives_2D(): + + ndim = 2 + + Px = sym.Symbol('Px') + Py = sym.Symbol('Py') + P0 = sym.Symbol('P0') + e33 = sym.Symbol('e33') + e31 = sym.Symbol('e31') + e15 = sym.Symbol('e15') + e00 = sym.Symbol('e00') + norm = sym.sqrt(Px**2 + Py**2) + P = np.array([Px, Py]) + I = np.eye(ndim) + p = P/norm + + # Piezoelectric tensor + eP = np.zeros((ndim,ndim,ndim),object) + for k,i,j in itertools.product(range(ndim), repeat=3): + eP[k,i,j] = (norm/P0)**3*(e33*p[i]*p[j]*p[k] + e31*(I[i,j]-p[i]*p[j])*p[k] + + e15*0.5*( (I[k,i]-p[k]*p[i])*p[j] + (I[k,j]-p[k]*p[j])*p[i] )) + + # derivative of eP wrt Px + deP_dPx = np.zeros_like(eP, object) + for k,i,j in itertools.product(range(ndim), repeat=3): + deP_dPx[k,i,j] = sym.diff(eP[k,i,j], Px) + + # derivative of eP wrt Py + deP_dPy = np.zeros_like(eP, object) + for k,i,j in itertools.product(range(ndim), repeat=3): + deP_dPy[k,i,j] = sym.diff(eP[k,i,j], Py) + + # calculate spon. strain + eps0_P = np.zeros((ndim,ndim), object) + for i,j in itertools.product(range(ndim), repeat=2): + eps0_P[i,j] = 3/2*e00*(norm/P0)**2*(p[i]*p[j] - 1/3*I[i,j]) + + # get its deivative wrt Px + deps0_P_dPx = np.zeros_like(eps0_P, object) + for i,j in itertools.product(range(ndim), repeat=2): + deps0_P_dPx[i,j] = sym.diff(eps0_P[i,j], Px) + + # get its deivative wrt Py + deps0_P_dPy = np.zeros_like(eps0_P, object) + for i,j in itertools.product(range(ndim), repeat=2): + deps0_P_dPy[i,j] = sym.diff(eps0_P[i,j], Py) + + + eP_val = sym.lambdify([Px, Py, P0, e33, e31, e15], eP, 'numpy') + deP_dPx_val = sym.lambdify([Px, Py, P0, e33, e31, e15], deP_dPx, 'numpy') + deP_dPy_val = sym.lambdify([Px, Py, P0, e33, e31, e15], deP_dPy, 'numpy') + + eps0_val = sym.lambdify([Px, Py, P0, e00], eps0_P, 'numpy') + deps0_P_dPx_val = sym.lambdify([Px, Py, P0, e00], deps0_P_dPx, 'numpy') + deps0_P_dPy_val = sym.lambdify([Px, Py, P0, e00], deps0_P_dPy, 'numpy') + + del eP, deP_dPx, deP_dPy, eps0_P, deps0_P_dPy, deps0_P_dPx + + return eP_val, deP_dPx_val, deP_dPy_val, eps0_val, deps0_P_dPx_val, deps0_P_dPy_val + +def Calculate_Derivatives_Land_Potential_2D(): + + Px = sym.Symbol('Px') + Py = sym.Symbol('Py') + P0 = sym.Symbol('P0') + + a0 = sym.Symbol('a0') + a1 = sym.Symbol('a1') + a2 = sym.Symbol('a2') + a3 = sym.Symbol('a3') + a4 = sym.Symbol('a4') + + Psi = a0 + a1/P0**2 * (Px**2 + Py**2) + a2/P0**4 * (Px**4 + Py**4) \ + + a3/P0**4 * Px**2 * Py**2 + a4/P0**6 * (Px**6 + Py**6) + + dPsi_dPx = sym.diff(Psi, Px) + dPsi_dPy = sym.diff(Psi, Py) + + Psi_val = sym.lambdify([Px, Py, P0, a0, a1, a2, a3, a4], Psi, 'numpy') + dPsi_dPx_val = sym.lambdify([Px, Py, P0, a1, a2, a3, a4], dPsi_dPx, 'numpy') + dPsi_dPy_val = sym.lambdify([Px, Py, P0, a1, a2, a3, a4], dPsi_dPy, 'numpy') + + del Psi, dPsi_dPx, dPsi_dPy + + return Psi_val, dPsi_dPx_val, dPsi_dPy_val diff --git a/fourier_frequencies.py b/fourier_frequencies.py new file mode 100644 index 0000000000000000000000000000000000000000..c4236715b39f2ab6d189adf3cc45405571ad9efe --- /dev/null +++ b/fourier_frequencies.py @@ -0,0 +1,39 @@ +""" +Calculates frequencies in Fourier space +Returns the vector of frequencies +""" + +import numpy as np + +def fourier_space(Nx, dx, Ny, dy, Nz=None, dz=None): + + if Nz == None: + + kx = 2.0*np.pi*np.fft.fftfreq(Nx, dx) + ky = 2.0*np.pi*np.fft.fftfreq(Ny, dy) + + kx_grid, ky_grid = np.meshgrid(kx,ky, indexing='ij') + + freq = np.array([kx_grid, ky_grid]) + + del kx, ky, kx_grid, ky_grid + + return freq + + else: + + kx = 2.0*np.pi*np.fft.fftfreq(Nx, dx) + ky = 2.0*np.pi*np.fft.fftfreq(Ny, dy) + kz = 2.0*np.pi*np.fft.fftfreq(Nz, dz) + + kx_grid, ky_grid, kz_grid = np.meshgrid(kx, ky, kz, indexing='ij') + + freq = np.array([kx_grid, ky_grid, kz_grid]) + + del kx, ky, kz, kx_grid, ky_grid, kz_grid + + return freq + + + + \ No newline at end of file diff --git a/green_operator.py b/green_operator.py new file mode 100644 index 0000000000000000000000000000000000000000..7ab9211c18b3230cc82e7c50584ce6219e4c3dc8 --- /dev/null +++ b/green_operator.py @@ -0,0 +1,45 @@ +""" +Calculates the Green's functions in Fourier space +""" + +import numpy as np +import itertools + +def Green_Operator_Piezoelectricity_2D(C0, e0, K0, freq): + + A_c = np.einsum('pijqx..., px..., qx... -> ijx...', C0, freq, freq) + A_e = np.einsum('piqx..., px..., qx... -> ix... ', e0, freq, freq) + A_d = np.einsum('pqx..., px..., qx... -> x... ', K0, freq, freq) + + A_d[0,0] = 1.0 + + invA_d = 1/A_d + + A = A_c + invA_d * np.einsum('ix..., jx... -> ijx... ', A_e, A_e) + + adjG = np.empty_like(A) + + adjG[0, 0] = A[1, 1] + adjG[1, 1] = A[0, 0] + adjG[0, 1] = -A[0, 1] + adjG[1, 0] = -A[1, 0] + + detG = A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0] + detG[0,0] = 1.0 + + invA = adjG/detG + + G_elas = np.zeros_like(C0) # Green op. for elasticity + for i,j,k,l in itertools.product(range(2), repeat=4): + G_elas[i,j,k,l] = 0.25*(invA[i,l]*freq[k]*freq[j]+invA[j,l]*freq[k]*freq[i]+invA[i,k]*freq[l]*freq[j]+invA[j,k]*freq[i]*freq[l]) + + G_piezo = np.zeros_like(e0) # Green op. for piezoelectricity + for k,i,j in itertools.product(range(2), repeat=3): + G_piezo[k,i,j] =( 1/2 * invA_d *freq[k]*(invA[i,0]*freq[j] + invA[j,0]*freq[i])*A_e[0] + + 1/2 * invA_d * freq[k]*(invA[i,1]*freq[j] + invA[j,1]*freq[i])*A_e[1] ) + + G_dielec = invA_d**2*(np.einsum('i...,ij...,j...',A_e,invA,A_e) - A_d)*np.einsum('ix..., jx... -> ijx... ', freq, freq) + + del A_c, A_e, A_d, invA_d, A, adjG, detG, invA + + return G_elas, G_piezo, G_dielec diff --git a/make_plot.py b/make_plot.py new file mode 100644 index 0000000000000000000000000000000000000000..ca6c1e0247937e8eafc5befeb43bf6abeea1dbc4 --- /dev/null +++ b/make_plot.py @@ -0,0 +1,67 @@ +""" +This script is used for postprocessing: + - plots polarization magnitude with its directions as an image. + - plots 6th-order Landau potential +""" + +import matplotlib.pyplot as plt +import numpy as np +from scipy import ndimage + +def Plot_Polarization(variable, tt, Nx, Ny, P1, P2, max_value): + + x,y = np.meshgrid(np.linspace(0,Nx,Nx),np.linspace(0,Ny,Ny)) + u = ndimage.rotate(P1, 90) + v = ndimage.rotate(P2, 90) + r2 = np.power(np.add(np.power(u,2),np.power(v,2)),0.5) + + widths = np.linspace(0, 1, x.size) + fig, ax = plt.subplots(figsize=(10, 10)) + + P_mag = np.sqrt(P1**2 + P2**2) + P_mag = ndimage.rotate(P_mag, 90) + #""" + im = ax.imshow(P_mag, cmap='jet', vmin=0, vmax=max_value) + ax.set_xlabel('x') + ax.set_ylabel('y') + + #""" + # colorbar + cax = fig.add_axes([ax.get_position().x1+0.01,ax.get_position().y0,0.02,ax.get_position().height]) + cb = plt.colorbar(im,cax=cax) + tick_font_size = 20 + cb.ax.tick_params(labelsize=tick_font_size) + + u = u/r2; v= v/r2 + + n = 14 + ax.quiver(x[::n,::n],y[::n,::n],u[::n,::n],v[::n,::n],linewidths=widths,scale=3.8, + scale_units='inches', color='white') + ax.axis('off') + fig.savefig(variable + str(tt), dpi=fig.dpi, bbox_inches='tight')#, pad_inches=0.5) + plt.close(fig) + + del u, v, P_mag + +def Plot_Potential(folder, P0, Nx, a0, a1, a2, a3, a4): + + from matplotlib import cm + + fig, ax = plt.subplots(subplot_kw={"projection": "3d"}) + + # Make data. + Px = np.linspace(-P0, P0, Nx) + Py = np.linspace(-P0, P0, Nx) + Px, Py = np.meshgrid(Px, Py) + Psi = a0 + a1/P0**2 * (Px**2 + Py**2) + a2/P0**4 * (Px**4 + Py**4) \ + + a3/P0**4 * (Px**2 * Py**2) + a4/P0**6 * (Px**6 + Py**6) + + + # Plot the surface. + ax.plot_surface(Px, Py, Psi, cmap=cm.coolwarm, + linewidth=0, antialiased=False) + + ax.zaxis.set_major_formatter('{x:.02f}') + + fig.savefig(folder + "/Landau_potential", dpi=fig.dpi, bbox_inches='tight') + plt.close(fig) \ No newline at end of file diff --git a/material_data.py b/material_data.py new file mode 100644 index 0000000000000000000000000000000000000000..8d83185ccb422ab7498a8f76cb835c26cef1a7e5 --- /dev/null +++ b/material_data.py @@ -0,0 +1,188 @@ +""" +Material data: + - store all material parameters here + - function Mat_Data_BTO_MD() contains material parameters for BTO + +Initial polarization can also be imported from here for different domain types: + - random -> polarization vector field starts with an initial + uniform distribution over [-0.1, +0.1] + - 90 -> polarization vector field starts with three 90° sharp interfaces + Note, here must be Nx=Ny + - 180 -> polarization vector field starts with a single 180° sharp interface +""" + +import numpy as np +import itertools + +# ---------------------------------------------------------------------------- +def full_3x3_to_Voigt_6_index_3D(i, j): + if i == j: + return i + return 6-i-j + +def full_3x3_to_Voigt_6_index_2D(i, j): + if i == j: + return i + return 3-i-j + + +def Voigt_to_full_Tensor_2D(C): + C = np.asarray(C) + C_out = np.zeros((2,2,2,2), dtype=float) + for i, j, k, l in itertools.product(range(2), range(2), range(2), range(2)): + Voigt_i = full_3x3_to_Voigt_6_index_2D(i, j) + Voigt_j = full_3x3_to_Voigt_6_index_2D(k, l) + C_out[i, j, k, l] = C[Voigt_i, Voigt_j] + return C_out + +def Voigt_to_full_Tensor_3D(C): + C = np.asarray(C) + C_out = np.zeros((3,3,3,3), dtype=float) + for i, j, k, l in itertools.product(range(3), range(3), range(3), range(3)): + Voigt_i = full_3x3_to_Voigt_6_index_3D(i, j) + Voigt_j = full_3x3_to_Voigt_6_index_3D(k, l) + C_out[i, j, k, l] = C[Voigt_i, Voigt_j] + return C_out +# ---------------------------------------------------------------------------- +def Mat_Data_BTO_D(): + + """ + https://reader.elsevier.com/reader/sd/pii/S0020768314000699?token=27FA350F30F884D46E6AA988FD02CD24C8E6CA0651D6635DA5E930EDCC7776CD8EA032BB95A0C94D524075D162599109&originRegion=eu-west-1&originCreation=20220329120452 + """ + + G = 1*12e-3 # J/m²=N/m Interfacial Energy + l_wall = 1.5E-9 # m: Thickness of interface: Length Scale + P0 = 0.26 # C/m²: Maximum polarization + mob = 26/75*1e3 # A/Vm: mobility beta^-1 + c_tet = 4.032 # angstrom + a_tet = 3.992 # angstrom + e00 = 2*(c_tet-a_tet)/(c_tet+2*a_tet) + + # ----------------------- MATERIAL PARAMETERS ---------------------------- + k_sep = 0.70 # Separation coefficient + k_grad = 0.35 # Gradient Energy coefficient + + # Parameters for Landau potential = k_sep*(G/l)*f(P) + a0 = 1 + a1 = -86/75 + a2 = -53/75 + a3 = 134/25 + a4 = 64/75 + l_coef = np.array([a0, a1, a2, a3, a4]) + + # piezoelectric tensor components + e31, e33, e15 = (-0.7, 6.7, 34.2) # e33 = 6.7 + e_piezo = np.array([e33, e31, e15]) + + # Elastic tensor comps + C11 = 22.2e10 + C12 = 11.1e10 + C44 = 6.1e10 + + C_cub = np.array([[C11, C12, 0],[C12, C11, 0],[0,0,C44]]) #*0.5 + + C0 = Voigt_to_full_Tensor_2D(C_cub) + + K11 = 19.5e-9 + K22 = 19.5e-9 + + K0 = np.array([[K11, 0],[0, K22]]) + + return G, l_wall, P0, mob, e00, k_sep, k_grad, l_coef, e_piezo, C0, C_cub, K0 + + +def Mat_Data_BTO_MD(): + + """ + https://reader.elsevier.com/reader/sd/pii/S0020768314000699?token=27FA350F30F884D46E6AA988FD02CD24C8E6CA0651D6635DA5E930EDCC7776CD8EA032BB95A0C94D524075D162599109&originRegion=eu-west-1&originCreation=20220329120452 + """ + + G = 12e-3 # J/m²=N/m Interfacial Energy + l_wall = 0.6E-9 # m: Thickness of interface: Length Scale + P0 = 0.26 # C/m²: Maximum polarization + mob = 2784#1857#26/75*1e3 # A/Vm: mobility beta^-1 + c_tet = 4.05 # angstrom + a_tet = 4 # angstrom + e00 = 2*(c_tet-a_tet)/(c_tet+2*a_tet) + + # ----------------------- MATERIAL PARAMETERS ---------------------------- + # Parameters for Landau potential = k_sep*(G/l)*f(P) + xi = 0.5 + phi_90 = 0.7 + + a0 = 1 + a1 = -(1-phi_90-2*xi**6)/(xi**2*(1-xi**4)) + a2 = -(-2*(1-phi_90)+3*xi**2+xi**6)/(xi**2*(1-xi**4)) + a3 = ((-2*xi**6+6*xi**4-3*xi**2+3*xi**2*phi_90+1-phi_90)/(xi**4*(1+xi**2))) #* 0.1 + a4 = (2*xi**2 + phi_90 - 1)/(xi**2*(1-xi**4)) + + l_coef = np.array([a0, a1, a2, a3, a4]) + + # Calibration coefficients + def integrand(x, a0, a1, a2, a4): + return np.sqrt(a0 + a1*x**2 + a2*x**4 + a4*x**6) + + from scipy.integrate import quad + Gamma, err = quad(integrand, -1, 1, args=(a0, a1, a2, a4)) + + k_grad = 0.5*(Gamma)**(-1) + k_sep = Gamma**(-1) + + # piezoelectric tensor components + e31, e33, e15 = (-0.7, 6.7, 34.2) # e33 = 6.7 + e_piezo = np.array([e33, e31, e15]) + + # Elastic tensor comps + C11 = 22.2e10 + C12 = 11.1e10 + C44 = 6.1e10 + + C_cub = np.array([[C11, C12, 0],[C12, C11, 0],[0,0,C44]]) #*0.1 + + C0 = Voigt_to_full_Tensor_2D(C_cub) + + K11 = 19.5e-9 + K22 = 19.5e-9 + + K0 = np.array([[K11, 0],[0, K22]]) + + return G, l_wall, P0, mob, e00, k_sep, k_grad, l_coef, e_piezo, C0, C_cub, K0 + + +def initial_polarization(Nx, Ny, domain_type): + + #np.random.seed(514) + + if domain_type == 'random': + Px = 0.1 * (2.0*np.random.rand(Nx, Ny) - 1.0) + Py = 0.1 * (2.0*np.random.rand(Nx, Ny) - 1.0) + + elif domain_type == '180': + Px = np.zeros((Nx,Ny)) + Py = np.ones((Nx,Ny)) + Py[int(Nx/2):,:]=-1 + + elif domain_type == '90': + + Px = np.zeros((Nx,Ny)) + Py = np.zeros((Nx,Ny)) + from scipy import linalg + first_row = np.zeros(Ny) + nn=Nx/2 + first_row[:int(nn)]=1 + first_row[2*int(nn)-1:3*int(nn)]=1 + + first_col = np.ones(Nx) + first_col[:int(nn)]=0 + first_col[2*int(nn)-1:3*int(nn)]=0 + + Py = linalg.toeplitz(first_col, first_row) + Px = (1-Py)*-1 + + Py = Py; Px=Px + + return np.array([Px, Py]) + + + + \ No newline at end of file diff --git a/solve_piezoelectricity.py b/solve_piezoelectricity.py new file mode 100644 index 0000000000000000000000000000000000000000..ffb4cf6dd5d251f741f0ea7deea1164f82034fc5 --- /dev/null +++ b/solve_piezoelectricity.py @@ -0,0 +1,79 @@ +""" +Interative Fourier spectral algorithm to solve the coupled piezoelectric constitutive equation. + +The method used here is based on the FFT-homogenization technique: + - introduces homogeneous C⁰(stiffness), e⁰(piezoelectric) and K⁰(dielectric) material tensors + these tensors are the volume average of heterogeous tensors C(P), e(P), K(p) + - Green's function for each tensor must be calculated first in Fourier space + - returns results for stresses, elastic strain, dielectric displacement and total electric field + +Variables used in Solve_Piezoelectricity_2D: + - C0, e0, K0 -> homogeneous material tensors, + - CP, eP, KP -> heterogeous material tensors (usually depent on the order parameter), + - G_elas, G_piezo, G_dielec -> Green's function for each tensor defined in Fourier space + - eps0 -> sponteneous polarization (depent on the order parameter), + - eps_ext -> applied strain, + - E_ext -> applied electric field, + - E_random -> random electric field associated with defects + - Nx, Ny -> grid points in each direction, + - P -> polarization (order parameter). +""" + +import numpy as np + +def Solve_Piezoelectricity_2D(C0,K0,e0,CP,KP,eP,G_elas,G_piezo,G_dielec,eps0,eps_ext,E_ext,E_random,Nx,Ny,P): + sig_norm = 1e-8 + D_norm = 1e-8 + # ------------------------------------------------------------------------ + eps_tot = np.zeros_like(eps0) + E = np.zeros_like(E_ext) + # ------------------------------------------------------------------------ + P_fft = np.fft.fft2(P,[Nx,Ny]) + eps0_fft = np.fft.fft2(eps0,[Nx,Ny]) + # ------------------------------------------------------------------------ + niter=100 + for itr in range(niter): + print(itr) + # -------------------------------------------------------------------- + eps_elas = eps_tot+eps_ext-eps0 + E_tot = E + E_ext + E_random + # -------------------------------------------------------------------- + sigma = np.einsum('ijklxy,klxy->ijxy',CP, eps_elas) - np.einsum('kijxy,kxy->ijxy',eP,E_tot) + D = np.einsum('ijkxy,jkxy->ixy',eP, eps_elas) + np.einsum('ijxy,jxy->ixy',KP, E_tot) + P + # convergence check ------------------------------------------------- + sig_norm_new = np.linalg.norm(sigma.ravel(), ord=2) + D_sum = np.sum(D, axis=0) + D_norm_new = np.linalg.norm(D_sum, ord=2) + err_s = np.abs((sig_norm_new - sig_norm) / sig_norm) + err_d = np.abs((D_norm_new - D_norm )/ D_norm) + print('err_s: ',err_s) + print('err_d: ',err_d) + if np.max(np.array([err_s,err_d])) < 1e-4: + break + D_norm = D_norm_new + sig_norm = sig_norm_new + # ------------------------------------------------------------------- + tau = sigma - np.einsum('ijklxy,klxy->ijxy',C0, eps_elas) + np.einsum('kijxy,kxy->ijxy',e0,E_tot) + rho = D - np.einsum('ijkxy,jkxy->ixy',e0, eps_elas) - np.einsum('ijxy,jxy->ixy',K0, E_tot) - P + # -------------------------------------------------------------------- + tau_fft = np.fft.fft2(tau,[Nx,Ny]) + rho_fft = np.fft.fft2(rho,[Nx,Ny]) + # -------------------------------------------------------------------- + alpha = tau_fft-np.einsum('ijklx...,klx...->ijx...',C0,eps0_fft) + beta = rho_fft + P_fft - np.einsum('ijkx...,jkx...->ix...',e0,eps0_fft) + # -------------------------------------------------------------------- + eps_tot_fft = -1*np.einsum('ijklx...,klx...->ijx...',G_elas,alpha) - np.einsum('kijx...,kx...->ijx...',G_piezo, beta) + eps_tot_fft[:,:,0,0]=0 + # -------------------------------------------------------------------- + E_fft = np.einsum('ijkx...,jkx...->ix...',G_piezo, alpha) + np.einsum('ijx...,jx...->ix...',G_dielec, beta) + E_fft[:,0,0]=0 + # -------------------------------------------------------------------- + eps_tot = np.fft.ifft2(eps_tot_fft,[Nx,Ny]).real + E = np.fft.ifft2(E_fft,[Nx,Ny]).real + # -------------------------------------------------------------------- + del eps_tot, E, P_fft, eps0_fft, tau, rho + del tau_fft, rho_fft, alpha, beta, E_fft, eps_tot_fft + # -------------------------------------------------------------------- + return sigma,D,eps_elas,E_tot + + diff --git a/store_hdf5.py b/store_hdf5.py new file mode 100644 index 0000000000000000000000000000000000000000..249f611e6a40bf3cb27f74a8969574eae0d22863 --- /dev/null +++ b/store_hdf5.py @@ -0,0 +1,86 @@ +""" +All data can be stored as .h5 file + +Write_to_HDF5 -> in this function you always need to give 3 NON KEYWORD ARGUMENTS: + + - folder -> directory to save results + - step -> the time step + - P -> polarization + +Write_to_HDF5 also contains KEYWORD ARGUMENTS + + - save data with the following keys and value: + + Stress=sigma + Elas_strain=eps_elas + Elec_disp=D + Elec_field=E + Elas_en=H_elas + Grad_en=H_grad + Sep_en=H_sep + Elec_en=H_elec + Piezo_en=H_piezo +""" + +import h5py + +def Write_to_HDF5(folder, step, P, **kwargs): + + hdf = h5py.File(folder + '/results.h5','a') + time = '/time_'+str(step) + + hdf.create_dataset('Polarization/Px'+str(time), data = P[0]) + hdf.create_dataset('Polarization/Py'+str(time), data = P[1]) + + for key, value in kwargs.items(): + + if key == 'Stress': + hdf.create_dataset('Stress/stress_XX'+str(time), data = value[0,0]) + hdf.create_dataset('Stress/stress_XY'+str(time), data = value[0,1]) + hdf.create_dataset('Stress/stress_YX'+str(time), data = value[1,0]) + hdf.create_dataset('Stress/stress_YY'+str(time), data = value[1,1]) + + elif key == 'Elas_strain': + hdf.create_dataset('Elastic strain/strain_XX'+str(time), data = value[0,0]) + hdf.create_dataset('Elastic strain/strain_XY'+str(time), data = value[0,1]) + hdf.create_dataset('Elastic strain/strain_YX'+str(time), data = value[1,0]) + hdf.create_dataset('Elastic strain/strain_YY'+str(time), data = value[1,1]) + + elif key == 'Elec_disp': + hdf.create_dataset('Diel. displacement/Dx'+str(time), data = value[0]) + hdf.create_dataset('Diel. displacement/Dy'+str(time), data = value[1]) + + elif key == 'Elec_field': + hdf.create_dataset('Electric field/Ex'+str(time), data = value[0]) + hdf.create_dataset('Electric field/Ey'+str(time), data = value[0]) + + elif key == 'Elas_en': + hdf.create_dataset('Elastic energy'+str(time), data = value) + + elif key == 'Grad_en': + hdf.create_dataset('Gradient energy'+str(time), data = value) + + elif key == 'Sep_en': + hdf.create_dataset('Landau energy'+str(time), data = value) + + elif key == 'Elec_en': + hdf.create_dataset('Electric energy'+str(time), data = value) + + elif key == 'Piezo_en': + hdf.create_dataset('Piezoelectric energy'+str(time), data = value) + + hdf.close() + + + + + + + + + + + + + +