1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine inversewavspd(xi,et,c,rho,a) 20! 21! Calculates the propagation wave speed in a material, up to its 21 22! constants. Subroutine for calcmatwavsps.f 23 24! Based on the procedure in: 25! C. Lane. The Development of a 2D Ultrasonic Array Inspection 26! for Single Crystal Turbine Blades. 27! Switzerland: Springer International Publishing, 2014. 28! 29! CARLO MONJARAZ TEC (CMT) 30! 31! INPUT: 32! 33! xi,et: values within a square domain between -1 and 1. 34! correlate in a unique way with 0<=phi<=pi and 35! -pi<=theta<=pi 36! 37! elas: c(3,3,3,3) - The elasticity vector 38! 39! rho: double - Density of the material 40! 41! 42! OUTPUT: 43! 44! a: inverse of the wave speed 45! 46 implicit none 47! 48 integer i,j,k,l,ier,matz,ndim 49! 50 real*8 c(3,3,3,3),rho,xn(3),cm(3,3,3),a,xi,et, 51 & cmm(3,3),dd,al(3),alz(3,3),fv1(3),fv2(3), 52 & theta,phi,pi,p3(3),v(3), 53 & speed 54! 55! 56! 57 pi=4.d0*datan(1.d0) 58! 59 theta=xi*pi 60 phi=(et+1.d0)*pi/2.d0 61! 62 do l=1,3 63 v(l)=0.d0 64 do k=1,3 65 cmm(k,l)=0.d0 66 do j=1,3 67 cm(j,l,k)=0.d0 68 enddo 69 enddo 70 enddo 71! 72 xn(1)=dcos(theta)*dsin(phi) 73 xn(2)=dsin(theta)*dsin(phi) 74 xn(3)=dcos(phi) 75! 76! c ------------ PER EAcH DIREcTION find wave speed----------------------- 77! 78 do l=1,3 79 do k=1,3 80 do i=1,3 81 do j=1,3 82 cm(l,k,i)=cm(l,k,i)+c(l,k,j,i)*xn(j) 83 enddo 84 enddo 85 enddo 86 enddo 87! 88 do k=1,3 89 do i=1,3 90 do l=1,3 91 cmm(k,i)=cmm(k,i)+cm(l,k,i)*xn(l) 92 enddo 93 enddo 94 enddo 95! 96 ndim=3 97 matz=1 98 ier=0 99! 100! ---------reset vars for EIGvALUES 101! 102 do j=1,3 103 al(j)=0.d0 104 fv1(j)=0.d0 105 fv2(j)=0.d0 106 do i=1,3 107 alz(j,i)=0.d0 108 enddo 109 enddo 110! 111 call rs(ndim,ndim,cmm,al,matz,alz,fv1,fv2,ier) 112! 113! ------normalizing eigenvectors to P vectors---------- 114! 115 dd=dsqrt(alz(1,3)**2+alz(2,3)**2+alz(3,3)**2) 116 p3(1)=alz(1,3)/dd 117 p3(2)=alz(2,3)/dd 118 p3(3)=alz(3,3)/dd 119! 120 do l=1,3 121 do k=1,3 122 cmm(k,l)=0.d0 123 do j=1,3 124 cm(j,l,k)=0.d0 125 enddo 126 enddo 127 enddo 128! 129 do l=1,3 130 do j=1,3 131 do i=1,3 132 do k=1,3 133 cm(l,j,i)=cm(l,j,i)+c(l,k,j,i)*xn(k); 134 enddo 135 enddo 136 enddo 137 enddo 138! 139 do l=1,3 140 do j=1,3 141 do i=1,3 142 cmm(l,j)=cmm(l,j)+cm(l,j,i)*p3(i); 143 enddo 144 enddo 145 enddo 146! 147 do j=1,3 148 do l=1,3 149 v(j)=v(j)+cmm(l,j)*p3(l) 150 enddo 151 enddo 152! 153 dd=dsqrt(v(1)**2+v(2)**2+v(3)**2) 154 speed=dd/dsqrt(rho*al(3)) 155c speed=dsqrt(dd/rho) 156! 157! inverse wave speed 158! 159 a=1.d0/speed 160! 161 return 162 end 163 164