晶体的单晶声速计算公式
Christoffel’s equation relates acoustic velocitics to crystallo-graphic direction and elastic constants.
$$
det|C_{ijkl}n_{l}n_{j}-\rho v^{2}\delta_{ik}|=0 \
det|F_{ik}n_{l}n_{j}-\rho v^{2}\delta_{ik}|=0 \
i,j,k,l=1,2,3; \
when\ i=k; then\ \delta=1; others\ \delta=0
$$
So, the Kelvin-Christoffel stiffness $\Gamma_{ik}$ is:
$$
\Gamma_{ik}=\sum_{j,l=1}^{n=3} C_{ijkl}n_{j}n_{l},\ i,k=1,2,3 \
\Gamma = \begin{bmatrix} \Gamma_{11} & \Gamma_{12} & \Gamma_{13} \ \Gamma_{12} & \Gamma_{22} & \Gamma_{23} \ \Gamma_{13} & \Gamma_{23} & \Gamma_{33}\end{bmatrix}
$$
For $n_{j}(\phi,\theta,\chi),n_{l}(\phi,\theta,\chi)$ :
$$
n_1=cos(\phi)cos(\theta)cos(\chi)-sin(\phi)sin(\chi) \
n_2=-cos(\phi)cos(\theta)sin(\chi)-sin(\phi)cos(\chi) \
n_3=-cos(\phi)sin(\theta)
$$
then, we can calculate the eigenvalues of matrix equation ($\lambda_i,i=1,2,3$).
$$
\lambda_1=\rho v_l^2 \
\lambda_2=\rho v_{s1}^2 \
\lambda_3=\rho v_{s2}^2
$$
Reference:
Sound velocity and elasticity of single-crystal forsterite to 16 GPa; DOI: 10.1029/96JB01266
The fortran code:
module singal_sound_velocity
contains
subroutine sing_sound_velocity(Celas, theta, phi, chi, V_sing)
! Celas: elastic constants
implicit none
real(kind=8) :: theta, phi ,chi, Celas(6,6)
real(kind=8) :: V_sing(3,3)
real(kind=8) :: n1,n2,n3
real(kind=8) :: f(3,3), Fv(3,3),E(3,3), namda(3,3)
integer :: n=3
real(kind=8) :: err
integer :: i,j
n1=cosd(theta)*cosd(phi)*cosd(chi)-sind(phi)*sind(chi)
n2=-cosd(theta)*cosd(phi)*sind(chi)-sind(phi)*cosd(chi)
n3=sind(theta)*cosd(phi)
f(1,1)=n1*n1*Celas(1,1)+n1*n2*Celas(1,6)+n1*n3*Celas(1,5) &
+n2*n1*Celas(6,1)+n2*n2*Celas(6,6)+n2*n3*Celas(6,5) &
+n3*n1*Celas(5,1)+n3*n2*Celas(5,6)+n3*n3*Celas(5,5)
f(2,2)=n1*n1*Celas(6,6)+n1*n2*Celas(6,2)+n1*n3*Celas(6,4) &
+n2*n1*Celas(2,6)+n2*n2*Celas(2,2)+n2*n3*Celas(2,4) &
+n3*n1*Celas(4,6)+n3*n2*Celas(4,2)+n3*n3*Celas(4,4)
f(3,3)=n1*n1*Celas(5,5)+n1*n2*Celas(5,4)+n1*n3*Celas(5,3) &
+n2*n1*Celas(4,5)+n2*n2*Celas(4,4)+n2*n3*Celas(4,3) &
+n3*n1*Celas(3,5)+n3*n2*Celas(3,4)+n3*n3*Celas(3,3)
f(1,2)=n1*n1*Celas(1,6)+n1*n2*Celas(1,2)+n1*n3*Celas(1,4) &
+n2*n1*Celas(6,6)+n2*n2*Celas(6,2)+n2*n3*Celas(6,4) &
+n3*n1*Celas(5,6)+n3*n2*Celas(5,2)+n3*n3*Celas(5,4)
f(1,3)=n1*n1*Celas(1,5)+n1*n2*Celas(1,4)+n1*n3*Celas(1,3) &
+n2*n1*Celas(6,5)+n2*n2*Celas(6,4)+n2*n3*Celas(6,3) &
+n3*n1*Celas(5,5)+n3*n2*Celas(5,4)+n3*n3*Celas(5,3)
f(2,3)=n1*n1*Celas(6,5)+n1*n2*Celas(6,4)+n1*n3*Celas(6,3) &
+n2*n1*Celas(2,5)+n2*n2*Celas(2,4)+n2*n3*Celas(2,3) &
+n3*n1*Celas(4,5)+n3*n2*Celas(4,4)+n3*n3*Celas(4,3)
Fv(1,1)=f(1,1)
Fv(2,2)=f(2,2)
Fv(3,3)=f(3,3)
Fv(1,2)=f(1,2)
Fv(1,3)=f(1,3)
Fv(2,1)=f(1,2)
Fv(3,1)=f(1,3)
Fv(2,3)=f(2,3)
Fv(3,2)=f(2,3)
err=1.0d-8
call eigenvector_calc(Fv,3,namda,err)
V_sing=namda
open(unit=89,file="velocity.out")
do i=1,3
write(89,*) (namda(i,j),j=1,3)
enddo
do i=1,3
write(89,*) (Fv(i,j),j=1,3)
enddo
end subroutine sing_sound_velocity
subroutine eigenvector_calc(a,n,e,eps)
!==============================================================
! Compute all eigenvalues: real symmetric matrix a(n,n,)
! a*x = lambda*x
! method: the basic QR method
! Alex G. (January 2010)
!--------------------------------------------------------------
! input ...
! a(n,n) - array of coefficients for matrix A
! n - dimension
! eps - convergence tolerance
! output ...
! e(n) - eigenvalues
! iter - number of iterations to achieve the tolerance
! comments ...
! kmax - max number of allowed iterations
!==============================================================
implicit none
integer :: n ,iter
real(kind=8) :: a(n,n), e(n,n), eps
real(kind=8) :: q(n,n), r(n,n), w(n), an, Ajnorm, sum, e0,e1
integer :: k, i, j, m
integer, parameter::kmax=1000
! initialization
q = 0.0
r = 0.0
e0 = 0.0
do k=1,kmax ! iterations
! step 1: compute Q(n,n) and R(n,n)
! column 1
an = Ajnorm(a,n,1)
r(1,1) = an
do i=1,n
q(i,1) = a(i,1)/an
end do
! columns 2,...,n
do j=2,n
w = 0.0
do m=1,j-1
! product q^T*a result = scalar
sum = 0.0
do i=1,n
sum = sum + q(i,m)*a(i,j)
end do
r(m,j) = sum
! product (q^T*a)*q result = vector w(n)
do i=1,n
w(i) = w(i) + sum*q(i,m)
end do
end do
! new a'(j)
do i =1,n
a(i,j) = a(i,j) - w(i)
end do
! evaluate the norm for a'(j)
an = Ajnorm(a,n,j)
r(j,j) = an
! vector q(j)
do i=1,n
q(i,j) = a(i,j)/an
end do
end do
! step 2: compute A=R(n,n)*Q(n,n)
a = matmul(r,q)
! egenvalues and the average eigenvale
sum = 0.0
do i=1,n
e(i,i) = a(i,i)
sum = sum+e(i,i)*e(i,i)
end do
e1 = sqrt(sum)
! print here eigenvalues
! write (*,201) (e(i),i=1,n)
!201 format (6f12.6)
! check for convergence
if (abs(e1-e0) < eps) exit
! prepare for the next iteration
e0 = e1
end do
iter = k
if(k == kmax) write (*,*)'The eigenvlue failed to converge'
end subroutine eigenvector_calc
end module
function Ajnorm(a,n,j)
implicit none
integer n, j, i
real(kind=8) :: a(n,n), Ajnorm
real(kind=8) :: sum
sum = 0.0
do i=1,n
sum = sum + a(i,j)*a(i,j)
end do
Ajnorm = sqrt(sum)
end function Ajnorm
