Doubt with reduction of mass and stiffness matrices in a system of 2 pdes

hello, I have this code:

clear all;

clc;

%---------Accelerators sensitivities---------------------------------------

sens0 = 1005; %Channel 0

sens1 = 962.54; %Channel 1

sens2 = 1005.11; %Channel 2

sens3 = 992.35; %Channel 3

fs=2048;

model=createpde(2);

%-----------Define dimensions of the plate-------------------------------

len = 1.22; % length in x direction

width =1.22; % Width in y direction

edges=[]; % Matrix to store the edges that will receive the boundary conditions

% Coordinates of the boundaries in x and y--------------------------------

coords_x = [0 len]';

coords_y = [0 width]';

sq=0.5; %Width of the inner square

% Create the geometry description matrix

gdm = [];

for k = 1 %:numel(coords_x)-1

for j = 1 %:numel(coords_y)-1

% Check if the coordinates are completely inside the external plate

%if coords_x(k) >= 0 && coords_x(k) = len && coords_y(j) >= 0 && coords_y(j) == width

% Add the vertical dimension

gdm = [gdm; 2, 12, coords_x(k), coords_x(k)+sq, coords_x(k+1)-sq, coords_x(k+1),coords_x(k+1),coords_x(k+1),coords_x(k+1), coords_x(k+1)-sq,coords_x(k)+sq, coords_x(k),coords_x(k),coords_x(k),coords_y(j), coords_y(j),coords_y(j),coords_y(j), coords_y(j)+sq, coords_y(j+1)-sq, coords_y(j+1),coords_y(j+1),coords_y(j+1),coords_y(j+1),coords_y(j+1)-sq,coords_y(j)+sq];

%end

end

end

%--------------------------------------------------------------------------

gdm=gdm'; % Format required by decsg function

%---------------Create the geometry--------------------------------------

[dl,bt] = decsg(gdm);

% Construct the geometry with the matrix dl so we keep the information of

% the edges

geometryFromEdges(model,dl);

% Generate the mesh and plot the geometry---------------------------------

msh = generateMesh(model, 'Hmax', 0.15); % Use 'Hmax' to control mesh density

figure %Visualise the model

pdemesh(model)

%---------Create the different regions of the model-----------------------

% Obtain nodes and elements of the mesh

[p,e,t] = meshToPet(model.Mesh);

%-----CREATE MODEL WITH THE GEOMETRIC DESCRIPTION IN DL--------------------------------------

% Definition of the constants

E = 3E9;

h_thick = 0.05;

nu = 0.3;

mass = 45; %kg/m2

D = E*h_thick^3/(12*(1-nu)^2);

% Now we create the PDE systems as symbolic equations_____________________

syms pres

syms u1(x,y,t) u2(x,y,t)

pdeeq = [-laplacian(u1,[x y])+u2==0; D*laplacian(u2,[x y])+ mass*diff(u1,t,t)==pres];

symcoeffs = pdeCoefficients(pdeeq,[u1,u2],'Symbolic',true);

c2=symcoeffs.c;

m2=symcoeffs.m;

%------Display the symbolic coefficients-----------------------------------

structfun(@disp,symcoeffs);

% Substitute a pression value in pres using the sub function. Since output

% of subs are symbolic objects, use pdecoefficienttodoubleto convert these

% coefficients to the double data type so they are valid for the pde toolbox

symcoeffs=subs(symcoeffs,pres,0);

symcoeffs=subs(symcoeffs);

coeffs=pdeCoefficientsToDouble(symcoeffs);

%------pass these coefficients to the pde model----------------------------

f = [0 0]'; %force the Force terms in both equations to be 0

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d,'c',...

coeffs.c,'a',coeffs.a,'f',f);

% INITIAL CONDITIONS ------------------------------------------------------

setInitialConditions(model,[0;0],[0;0]);

ac=DATA;

% -------Create vector time, not starting in 0

num_samples = size(DATA, 1);

dt = 1 / fs; % interval between measurements

t = (1:num_samples) * dt; % vector of times

% TRANSFORM ACCELERATIONS INTO DISPLACEMENTS_______________________________

acvedi=AccVelDis(sens0,sens1,sens2,ac,fs,t,sens3);

dis=acvedi{3};

t=[0,t]; % I start the time at 0

% Boundary conditions in the faces using the displacement vectors

in=0; %counter variable

Edges=[10 4 6 12]; %Vector with faces numbers

% rest of the edges

allEdges = 1:model.Geometry.NumEdges;

neumannEdges = [1 2 3 5 7 8 9 11];

for i = 1:numel(Edges)

in=in+1;

edge = Edges(i);

dis1=dis(:,in);

applyBoundaryCondition(model, "neumann", "Edge", edge, ...

"g", @(location,state) [mybc1(location,state,dis1,t); 0], ...

"q", [0 0; 1E6 0]);

end

tim=[2:99];

res=solvepde(model, tim);

% Access the solution at the nodal locations

sol_ver_spring06=res.NodalSolution;

when I use FEM = assembleFEMatrices(model,'KMFA') I obtain a k matrix that is nodes*2 x nodes*2 and likewise with the M matrix. I guess that they correspond to my variables U1 and U2 but I guess I cannot use those matrices to calculate the natural frequencies of the system? I have read about using the Schur complement to eliminate u2 from the block matrix system and to rewrite the system in a way that only involves u1. When I use spy (FEM.K) I have values in the 1st and 4th cuadrant of the matrix, whereas if I use spy(FEM.M) the only one with values is the 3rd cuadrant (it makes sense as it just affects the second variable u2):

M11 = FEM.M(1:nodes, 1:nodes);

M12 = FEM.M(1:nodes, nodes+1:end);

M21 = FEM.M(nodes+1:end, 1:nodes);

M22 = FEM.M(nodes+1:end, nodes+1:end);

K11 = FEM.K(1:nodes, 1:nodes);

K12 = FEM.K(1:nodes, nodes+1:end);

K21 =FEM.K(nodes+1:end, 1:nodes);

K22 = FEM.K(nodes+1:end, nodes+1:end);

% Schur complement to reduce the system

M_r = M11 - M12 * (M22 \ M21);

K_r = K11 - K12 * (K22 \ K21);

K_r seems to be ok as now the dimensions are nodesxnodes and M_r has the right dimensions but as just the 3rd cuadrant has values, M_r is empty.

Can I follow any other approach?

Thanks very much for your help

Answer 1

[UDAYA PEDDIRAJU]

Hi Jorge,

To address the issue with the mass matrix ( M_r ) being empty, you might want to check the following:

1. Verify that ( M_{12} ) and ( M_{21} ) are not zero matrices. If they are, the Schur complement will yield an empty matrix.
2. Ensure that the boundary conditions applied do not lead to a trivial solution, which could result in zero entries in the mass matrix.
3. If the Schur complement method does not yield the desired results, consider using modal analysis directly on the full system or applying a different reduction technique, such as Guyan reduction.

I am attaching my code here as an example:

% Check if M12 and M21 are non-zero
if all(M12(:) == 0) || all(M21(:) == 0)
    error('M12 or M21 is a zero matrix, check your boundary conditions.');
end

% Calculate the Schur complement
M_r = M11 - M12 * (M22 \ M21);
K_r = K11 - K12 * (K22 \ K21);

% Display the results
disp('Reduced Mass Matrix M_r:');
disp(M_r);
disp('Reduced Stiffness Matrix K_r:');
disp(K_r);

[Jorge Garcia Garcia]

also the problem does not seem to be well constrained, but I am not sure any longer how to achieve it. I need the edges to move following the same displacement that a displacement vector, different in every edge, therefore I cannot apply Dirichlet as it says when I run is not a Dirichlet boundary conditions. Any suggestion would be appreciated

Answer 2

[Jorge Garcia Garcia]

it makes sense that M12 and M11 are empty as there is just a m component in the second partial diferential equation. I will check the Guyan reduction.