Thin CCFF Plate bending

Question

abdallah samaha 2019-11-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/490856-thin-ccff-plate-bending

Hi, I have got this code to measure the FRF and mode shapes of thin plate under a point force with a boundry condition of CCSS(clamped on two sides and simply suppourted on the other two). I need to change it to the boundry conditions of CCFF(clamped on two sides and free on the other two). I also need to change tho point force to an area force. Can someone please help with this.

clc
clear all
close all
%% General Parameters:
%==========================================================================
% Frequency Range:
df0 = 1;
f0 = 1:df0:2000;
% Material Parameters:
E0 = 5.5e8;
rho0 = 7800;
h0 = 0.003;
nu0 = 0.33;
eta0 = 0.01;
% Dimensions in m:
L0 = [0.8,0.4];
% Number of Modes:
nmodes = 16;
% Point Force and Response Location in m:
F0 = 1;
xf0 = [0.6,0.3];
xr0 = [0.6,0.3];
%==========================================================================
%% Solving Characteristic Equation:
%==========================================================================
fun0 = @(beta) cos(beta)*cosh(beta)-1;
fund0 = @(beta) -sin(beta)*cosh(beta)+cos(beta)*sinh(beta);
fun1 = @(beta) beta-fun0(beta)./fund0(beta);
beta0 = ((1:4)+1/2)*pi;
beta1 = zeros(size(beta0));
for j=1:4
    beta00 = beta0(j);
    for k=1:1000
        beta1(j) = fun1(beta00);
        if(abs(fun0(beta1(j)))<1e-10)
            break
        else
            beta00 = beta1(j);
        end
    end
end
if(nmodes > 4)
    beta1 = [beta1,((5:nmodes)+1/2)*pi];
end
%==========================================================================
%% Vibrational Mode Shapes:
%==========================================================================
om0 = 2*pi*f0;
% Clamped Boundary Condition in x-direction:
k1 = beta1/L0(1);
gamma1 = @(m) (cos(k1(m)*L0(1))-cosh(k1(m)*L0(1)))./(sin(k1(m)*L0(1))-sinh(k1(m)*L0(1)));
mode1 = @(x,m) (cos(k1(m)*x)-gamma1(m)*sin(k1(m)*x))-cosh(k1(m)*x).*(1-gamma1(m)*tanh(k1(m)*x));
mode2 = @(x,m) (cos(k1(m)*x)-sin(k1(m)*x))+(sin(k1(m)*L0(1))-cos(k1(m)*L0(1))).*exp(-k1(m)*(L0(1)-x))-exp(-k1(m)*x);
modedd1 = @(x,m) k1(m)^2*(-(cos(k1(m)*x)-gamma1(m)*sin(k1(m)*x))-cosh(k1(m)*x).*(1-gamma1(m)*tanh(k1(m)*x)));
modedd2 = @(x,m) k1(m)^2*(-(cos(k1(m)*x)-sin(k1(m)*x))+(sin(k1(m)*L0(1))-cos(k1(m)*L0(1))).*exp(-k1(m)*(L0(1)-x))-exp(-k1(m)*x));
% Simply Supported Boundary Condition in y-direction:
k2 = (1:nmodes)*pi/L0(2);
mode3 = @(y,n) sin(k2(n)*y);
modedd3 = @(y,n) -k2(n)^2*sin(k2(n)*y);
%==========================================================================
%% Graphic Output:
%==========================================================================
x1 = linspace(0,L0(1),1000);
y1 = linspace(0,L0(2),1000);
% Clamped Vibrational Modes in x-direction:
figure(1)
subplot(1,2,1)
for k=1:4
    if(k<=5)
        plot(x1,-mode1(x1,k),'b-','LineWidth',2)
        hold on
    else
        plot(x1,-mode2(x1,k),'b-','LineWidth',2)
        hold on
    end
end
hold off
ylim([-2,2])
grid('on')
% Simply Supported Vibrational Modes in y-direction:
subplot(1,2,2)
for k=1:4
    plot(y1,mode3(y1,k),'b-','LineWidth',2)
    hold on
end
hold off
ylim([-2,2])
grid('on')
%==========================================================================
 
%% Stiffness and Mass:
%==========================================================================
D0 = E0*(1+1j*eta0)*h0^3/(12*(1-nu0^2));
m0 = rho0*h0;
% Initialize all possible modes:
[I0,J0] = meshgrid(1:nmodes,1:nmodes);
I0 = I0(:);
J0 = J0(:);
% Mass Matrix
M0 = m0*(prod(L0)/2)*speye(nmodes^2);
 
% Stiffness Matrix
K01 = diag(D0*(k1(I0).^4+k2(J0).^4)*(prod(L0)/2));
 
K02x = zeros(nmodes^2,1);
K02y = zeros(nmodes^2,1);
x0 = linspace(0,L0(1),1000);
y0 = linspace(0,L0(2),1000);
% clamped boundary condition:
k01 = @(x,m,p) modedd1(x,m).*mode1(x,p);   % m < 5 & p < 5
k02 = @(x,m,p) modedd1(x,m).*mode2(x,p);   % m < 5 & p >= 5
k03 = @(x,m,p) modedd2(x,m).*mode1(x,p);   % m >= 5 & p < 5
k04 = @(x,m,p) modedd2(x,m).*mode2(x,p);   % m >= 5 & p >= 5
% simply supported boundary condition:
k05 = @(y,n,q) modedd3(y,n).*mode3(y,q);
% clamped boundary condition in x-direction:
count0 = 1;
for k=1:nmodes
    for l=1:nmodes
        if(k<5 && l<5)
            K02x(count0) = trapz(x0,k01(x0,k,l));
            count0 = count0+1;
        elseif(k<5 && l>=5)
            K02x(count0) = trapz(x0,k02(x0,k,l));
            count0 = count0+1;
        elseif(k>=5 && l<5)
            K02x(count0) = trapz(x0,k03(x0,k,l));
            count0 = count0+1;
        elseif(k>=5 && l>=5)
            K02x(count0) = trapz(x0,k04(x0,k,l));
            count0 = count0+1;
        end
    end
end
% simply supported boundary condition in y-direction:
count0 = 1;
for k=1:nmodes
    for l=1:nmodes
        K02y(count0) = trapz(y0,k05(y0,k,l));
        count0 = count0+1;
    end
end
% Fusion of x- and y-direction:
K02x = reshape(K02x,nmodes,nmodes);
K02y = reshape(K02y,nmodes,nmodes);
K02x = reshape(K02x,nmodes,nmodes);
K02x = kron(K02x,ones(nmodes));
K02y = reshape(K02y,nmodes,nmodes);
K02y = repmat(K02y,nmodes,nmodes);
% Total stiffness matrix:
K00 = sparse(K01+2*D0*sparse(K02x.*K02y));
K0 = (K00+K00.')/2;
%==========================================================================
%% Point Force and Response:
%==========================================================================
load0 = zeros(nmodes^2,1);
resp0 = zeros(nmodes^2,1);
count0 = 1;
for k=1:numel(I0)
    if(I0(k)<5)
        load0(count0) = F0*mode1(xf0(1),I0(k)).*mode3(xf0(2),J0(k));
        resp0(count0) = mode1(xr0(1),I0(k)).*mode3(xr0(2),J0(k));
        count0 = count0+1;
    else
        load0(count0) = F0*mode2(xf0(1),I0(k)).*mode3(xf0(2),J0(k));
        resp0(count0) = mode2(xr0(1),I0(k)).*mode3(xr0(2),J0(k));
        count0 = count0+1;
    end
end
% Solution of equation system
w0 = zeros(size(f0));
for k=1:numel(f0)
    Mat0 = K0-om0(k)^2*M0;
    C0 = Mat0\load0;
    w0(k) = 1j*om0(k)*sum(C0.*resp0);
end
%==========================================================================
%% Graphic Output:
%==========================================================================
figure(1)
subplot(2,1,1)
plot(f0,10*log10((abs(w0)/5e-8).^2/2),'b-','LineWidth',2)
grid('on')
subplot(2,1,2)
plot(f0,(180/pi)*angle(w0),'b-','LineWidth',2)
grid('on')
%==========================================================================