Dec-26-2021, 09:04 AM
I am trying to find the solution to the equation. Det(matrix A)=0, var:omega
###################################################################
###################################################################
from sympy import Function, symbols,atan, sin, cos, pi, Matrix, solve
A11, A12, A22, A66, B11, B12, B22, B66, C11, C12, C22, C66, D11, D12, D22, D66, E11, E12, E22, E66, G11, G12, G22, G66, F44, F55, I1, I2, I3, I4, I5, I6\
=symbols('A11 A12 A22 A66 B11 B12 B22 B66 C11 C12 C22 C66 D11 D12 D22 D66 E11 E12 E22 E66 G11 G12 G22 G66 F44 F55 I1 I2 I3 I4 I5 I6')
omega =symbols('omega')
# from sympy.abc import x, y, t
u, U, v, V, w, W ,phix, Phix, phiy, Phiy, E, rho, muy, fz, gz= symbols('u U v V w W phix Phix phiy Phiy E rho muy fz gz', cls=Function)
from scipy.integrate import quad
###############################################################################
h=0.002
a=20*h
b=a
#mode
m=1
n=1
k0=0
k1=0
###############################################################################
def fz(z):
return h*(atan(3*z/h)+sin(2*pi*z/h)/3/pi)/5
def gz(z):
return (1/(1+2*(z**3)/(h**3))+cos(3*pi*z/h))/5
###############################################################################
Em=70*(10**9)
Ec=380*(10**9)
rhom=2702
rhoc=3800
muym=0.3
muyc=0.3
alpha=0
p=1
###############################################################################
def E(z):
return (Ec-Em)*((0.5+z/h)**p)+Em-alpha*(Ec+Em)/2
def rho(z):
return (rhoc-rhom)*((0.5+z/h)**p)+rhom-alpha*(rhoc+rhom)/2
def muy(z):
return (muyc-muym)*((0.5+z/h)**p)+muym-alpha*(muyc+muym)/2
###############################################################################
#Ma tran do cung Q
def Q11(z):
return E(z)/(1-muy(z)**2)
def Q12(z):
return E(z)*muy(z)/(1-muy(z)**2)
def Q22(z):
return E(z)/(1-muy(z)**2)
def Q66(z):
return E(z)/(1+muy(z))/2
###############################################################################
A11= quad(lambda z: Q11(z), -h/2, h/2)[0]
A12= quad(lambda z: Q12(z), -h/2, h/2)[0]
A22= quad(lambda z: Q22(z), -h/2, h/2)[0]
A66= quad(lambda z: Q66(z), -h/2, h/2)[0]
B11= quad(lambda z: Q11(z)*z, -h/2, h/2)[0]
B12= quad(lambda z: Q12(z)*z, -h/2, h/2)[0]
B22= quad(lambda z: Q22(z)*z, -h/2, h/2)[0]
B66= quad(lambda z: Q66(z)*z, -h/2, h/2)[0]
C11= quad(lambda z: Q11(z)*fz(z), -h/2, h/2)[0]
C12= quad(lambda z: Q12(z)*fz(z), -h/2, h/2)[0]
C22= quad(lambda z: Q22(z)*fz(z), -h/2, h/2)[0]
C66= quad(lambda z: Q66(z)*fz(z), -h/2, h/2)[0]
D11= quad(lambda z: Q11(z)*z*z, -h/2, h/2)[0]
D12= quad(lambda z: Q12(z)*z*z, -h/2, h/2)[0]
D22= quad(lambda z: Q22(z)*z*z, -h/2, h/2)[0]
D66= quad(lambda z: Q66(z)*z*z, -h/2, h/2)[0]
E11= quad(lambda z: Q11(z)*fz(z)*z, -h/2, h/2)[0]
E12= quad(lambda z: Q12(z)*fz(z)*z, -h/2, h/2)[0]
E22= quad(lambda z: Q22(z)*fz(z)*z, -h/2, h/2)[0]
E66= quad(lambda z: Q66(z)*fz(z)*z, -h/2, h/2)[0]
G11= quad(lambda z: Q11(z)*fz(z)*fz(z), -h/2, h/2)[0]
G12= quad(lambda z: Q12(z)*fz(z)*fz(z), -h/2, h/2)[0]
G22= quad(lambda z: Q22(z)*fz(z)*fz(z), -h/2, h/2)[0]
G66= quad(lambda z: Q66(z)*fz(z)*fz(z), -h/2, h/2)[0]
F44= quad(lambda z: Q66(z)*gz(z)*gz(z), -h/2, h/2)[0]
F55= quad(lambda z: Q66(z)*gz(z)*gz(z), -h/2, h/2)[0]
I1=quad(lambda z: rho(z), -h/2, h/2)[0]
I2=quad(lambda z: rho(z)*z, -h/2, h/2)[0]
I3=quad(lambda z: rho(z)*z*z, -h/2, h/2)[0]
I4=quad(lambda z: rho(z)*fz(z), -h/2, h/2)[0]
I5=quad(lambda z: rho(z)*fz(z)*z, -h/2, h/2)[0]
I6=quad(lambda z: rho(z)*fz(z)*fz(z), -h/2, h/2)[0]
a11=(a*A11*b**2*m**2*pi**2+a**3*A66*n**2*pi**2)/(a**3*b**2);
a12=(a**2*A12*b*m*n*pi**2+a**2*A66*b*m*n*pi**2)/(a**3*b**2);
a13=((-b**2)*B11*m**3*pi**3-a**2*B12*m*n**2*pi**3-2*a**2*B66*m*n**2*pi**3)/(a**3*b**2);
a14=((b**2*C11*m**2+a**2*C66*n**2)*pi**2)/(a**2*b**2);
a15=((C12+C66)*m*n*pi**2)/(a*b);
a21=(a*A12*b**2*m*n*pi**2+a*A66*b**2*m*n*pi**2)/(a**2*b**3);
a22=(A66*b**3*m**2*pi**2+a**2*A22*b*n**2*pi**2)/(a**2*b**3);
a23=((-b**2)*B12*m**2*n*pi**3-2*b**2*B66*m**2*n*pi**3-a**2*B22*n**3*pi**3)/(a**2*b**3);
a24=((C12+C66)*m*n*pi**2)/(a*b);
a25=((b**2*C66*m**2+a**2*C22*n**2)*pi**2)/(a**2*b**2);
a31=((-a)*b**4*B11*m**3*pi**3-a**3*b**2*B12*m*n**2*pi**3-2*a**3*b**2*B66*m*n**2*pi**3)/(a**4*b**4);
a32=((-a**2)*b**3*B12*m**2*n*pi**3-2*a**2*b**3*B66*m**2*n*pi**3-a**4*b*B22*n**3*pi**3)/(a**4*b**4);
a33=(1/(a**4*b**4))*(a**4*b**4*k0+a**2*b**4*k1*m**2*pi**2+a**4*b**2*k1*n**2*pi**2+b**4*D11*m**4*pi**4+2*a**2*b**2*D12*m**2*n**2*pi**4+4*a**2*b**2*D66*m**2*n**2*pi**4+a**4*D22*n**4*pi**4);
a34=-((m*(b**2*E11*m**2+a**2*(E12+2*E66)*n**2)*pi**3)/(a**3*b**2));
a35=-((n*(b**2*(E12+2*E66)*m**2+a**2*E22*n**2)*pi**3)/(a**2*b**3));
a41=(a*b**2*C11*m**2*pi**2+a**3*C66*n**2*pi**2)/(a**3*b**2);
a42=(a**2*b*C12*m*n*pi**2+a**2*b*C66*m*n*pi**2)/(a**3*b**2);
a43=((-b**2)*E11*m**3*pi**3-a**2*E12*m*n**2*pi**3-2*a**2*E66*m*n**2*pi**3)/(a**3*b**2);
a44=(a*b**2*G11*m**2*pi**2+a**3*(b**2*F55+G66*n**2*pi**2))/(a**3*b**2);
a45=((G12+G66)*m*n*pi**2)/(a*b);
a51=(a*b**2*C12*m*n*pi**2+a*b**2*C66*m*n*pi**2)/(a**2*b**3);
a52=(b**3*C66*m**2*pi**2+a**2*b*C22*n**2*pi**2)/(a**2*b**3);
a53=((-b**2)*E12*m**2*n*pi**3-2*b**2*E66*m**2*n*pi**3-a**2*E22*n**3*pi**3)/(a**2*b**3);
a54=((G12+G66)*m*n*pi**2)/(a*b);
a55=(b**3*G66*m**2*pi**2+a**2*(b**3*F44+b*G22*n**2*pi**2))/(a**2*b**3);
A=Matrix([[a11-I1*omega,a12,a13+(I2*m*pi*omega)/a,a14-I4*omega,a15],\
[a21,a22-I1*omega,a23+(I2*n*pi*omega)/b,a24,a25-I4*omega],\
[a31+(I2*m*pi*omega)/a,a32+(I2*n*pi*omega)/b,a33-((a**4*b**4*I1+a**2*b**4*I3*m**2*pi**2+a**4*b**2*I3*n**2*pi**2)*omega)/(a**4*b**4),a34+(I5*m*pi*omega)/a,a35+(I5*n*pi*omega)/b],\
[a41-I4*omega,a42,a43+(I5*m*pi*omega)/a,a44-I6*omega,a45],\
[a51,a52-I4*omega,a53+(I5*n*pi*omega)/b,a54,a55-I6*omega]])
DET=A.det()
ng=solve(DET, omega)
print(ng)