Commit df8f6612 authored by Andreas Zilian's avatar Andreas Zilian
Browse files

Added Stability analysis class.

parent c0f5691a
classdef FeAnalysisStability < mafe.FeAnalysis
% stability analysis of finite element problem
properties
X = []; % matrix of buckling modes, columns-wise
L = []; % vector of buckling load factors, eigenvalues
end
methods
%% constructor
function obj = FeAnalysisStability(fep)
if nargin > 0
obj.fep = fep;
end
end
%% build and compute the associated eigenproblem
function analyse(self)
% make sure the problem is initialised
self.fep.init();
% STEP 1: Perform linear analysis
% get system elastic stiffness matrix
K = self.fep.assembleSystemStiffnessMatrixElastic();
% get force vector
p = self.fep.assembleSystemForceVector();
% apply constraints
[f ] = self.fep.applyConstraintsNeumann(p);
[f, dbc, act] = self.fep.applyConstraintsDirichlet(K, f);
% allocate solution vector
u = zeros(self.fep.ndofs,1);
% solve system of equations (active indices only): elastic response, 1st order
u(act) = K(act,act)\f(act);
u(dbc) = f(dbc);
% compute reaction forces
r = K * u - p;
% feedback of solution values (u and r) to nodal dofs
self.fep.updateNodalDofs(u, r);
% postprocess elements
self.fep.postprocessElements(u);
% STEP 2: Perform linear stability analysis
% get system geometric stiffness matrix
G = self.fep.assembleSystemStiffnessMatrixGeometric();
% solve eigenvalue problem
% (a) compute eigenvalues lambda
sfac = 1./max( [ normest(K,1e-4), normest(G,1e-4) ] ); % scale just for polyeig
self.L = polyeig(K(act,act)*sfac, G(act,act)*sfac);
% (b) compute eigenvectors
self.X = zeros(self.fep.ndofs,length(self.L));
for ii = 1:length(self.L)
lambda = self.L(ii);
% check for invalid eigenvalues
if abs(lambda) == Inf % these are typically the membrane-only modes
continue
end
% create matrix
A = lambda * G(act,act) + K(act,act);
% compute eigenvector by singular value decomposition
% TODO: check if this is generalisable
[U,S,V] = svds(A, 1, 'smallest');
self.X(act,ii) = V(:,end);
% scale computed eigenvector
self.X(:,ii) = 1./sqrt(self.X(:,ii)'*K*self.X(:,ii)) * self.X(:,ii); % normalize mode wrt K
end
% (c) sort eigenvalues/eigenvectors in ascending order
[Lsorted, permutation] = sort( abs(self.L) );
self.X = self.X(:,permutation);
self.L = Lsorted;
end
%% put selected mode to dof solution
function putModeToDof(self, mode)
if nargin < 1
mode = 1;
end
lambda = self.L(mode);
u = 1 * self.X(:,mode);
u = 1./norm(u,2) * u;
% compute reaction forces
K = self.fep.assembleSystemStiffnessMatrix();
f = self.fep.assembleSystemForceVector();
r = K * u - f;
% feedback of solution values (u and r) to nodal dofs
self.fep.updateNodalDofs(u, r);
% postprocess elements
self.fep.postprocessElements(u);
end
end
end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment