Loading green_operator.py 0 → 100644 +45 −0 Original line number Diff line number Diff line """ 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 make_plot.py 0 → 100644 +67 −0 Original line number Diff line number Diff line """ 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 material_data.py 0 → 100644 +188 −0 Original line number Diff line number Diff line """ 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 solve_piezoelectricity.py 0 → 100644 +79 −0 Original line number Diff line number Diff line """ 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 store_hdf5.py 0 → 100644 +86 −0 Original line number Diff line number Diff line """ 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() Loading
green_operator.py 0 → 100644 +45 −0 Original line number Diff line number Diff line """ 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
make_plot.py 0 → 100644 +67 −0 Original line number Diff line number Diff line """ 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
material_data.py 0 → 100644 +188 −0 Original line number Diff line number Diff line """ 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
solve_piezoelectricity.py 0 → 100644 +79 −0 Original line number Diff line number Diff line """ 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
store_hdf5.py 0 → 100644 +86 −0 Original line number Diff line number Diff line """ 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()
