Commit 713c3ab7 authored by Andreas Zilian's avatar Andreas Zilian
Browse files

Added initial version of flat shell.

parent df8f6612
classdef FlatShell2D < mafe.Qu4Element2D
% Plate (Reissner theory) + Membrane finite element in the x-y plane (2D)
properties
membr;
plate;
idx_membr = [1:8];
idx_plate = [9:20];
end
methods
%% constructor
function obj = FlatShell2D(varargin)
obj = obj@mafe.Qu4Element2D(varargin{:});
obj.plate = mafe.Plate2D(varargin{:});
obj.membr = mafe.Membrane2D(varargin{:});
obj.sect.state = 'plane_stress';
obj.internal = zeros(8,4);
end
% ======================================================================
%% activation of dofs at the nodes of the element
function activateDofTypes(self)
ineed = [ mafe.DofType.Disp1 mafe.DofType.Disp2 mafe.DofType.Disp3 mafe.DofType.Rota1 mafe.DofType.Rota2 ];
for i = 1:length(self.node)
self.node(i).dof = [self.node(i).dof ineed];
end
end
%% provide element index vector of dof
function idx = getDofSysIndices(self)
disp1 = find( [self.node.dof] == mafe.DofType.Disp1 );
disp2 = find( [self.node.dof] == mafe.DofType.Disp2 );
disp3 = find( [self.node.dof] == mafe.DofType.Disp3 );
rota1 = find( [self.node.dof] == mafe.DofType.Rota1 );
rota2 = find( [self.node.dof] == mafe.DofType.Rota2 );
ndidx = [ self.node.idx ];
idx = [ ndidx(disp1), ndidx(disp2), ndidx(disp3), ndidx(rota1), ndidx(rota2) ];
end
% ======================================================================
%% calculation of element stiffness matrix
function ele_mat = stiffness(self)
% global element matrix
ele_mat = self.stiffnessElastic() + self.stiffnessGeometric();
end
% calculation of elastic element stiffness matrix
function ele_mat = stiffnessElastic(self)
Ke_membr = self.membr.stiffness();
Ke_plate = self.plate.stiffness();
ele_mat = zeros(20,20);
ele_mat(self.idx_membr,self.idx_membr) = Ke_membr;
ele_mat(self.idx_plate,self.idx_plate) = Ke_plate;
end
% calculation of geometric element stiffness matrix
function ele_mat = stiffnessGeometric(self)
% overall matrix
ele_mat = zeros(20,20);
% sub-matrix of the 2nd order term
rotandof = 2*4;
ele_mat_rr = zeros(rotandof,rotandof);
% Gauss integration parameters
[GC, GW] = mafe.GaussCoordWeights(2);
% normal (membrane) forces
M = [ mean(self.internal(1,:)), mean(self.internal(3,:));
mean(self.internal(3,:)), mean(self.internal(2,:)) ];
% numerical integration of weak form (end order term)
for ii = 1:length(GC) % looping over gauss points in z1 direction
for jj = 1:length(GC) % looping over gauss points in z2 direction
% gauss coordinates
z = [ GC(ii), GC(jj) ];
% gauss weights
w = GW(ii) * GW(jj);
% shape/ansatz functions
J = self.J(z); % Jacobi matrix
N = self.N(z); % local derivative
detJ = det(J); % Jacobian
% integration factor
factor = w * detJ * self.sect.t;
% strain approach (in global coordinates)
A = [ N, zeros(1,4) ;
zeros(1,4), N ];
% internal virtual work expressions
ele_mat_rr = ele_mat_rr + factor * A' * M * A;
end
end
ele_mat(13:20,13:20) = ele_mat_rr;
end
%% calculation of element mass matrix
function ele_mat = mass(self)
%
ele_mat = zeros(20,20);
%
ele_mat(self.idx_membr,self.idx_membr) = self.membr.mass();
ele_mat(self.idx_plate,self.idx_plate) = self.plate.mass();
end
%% calculation of element damping matrix
function ele_mat = damping(self)
%
ele_mat = zeros(20,20);
%
ele_mat(self.idx_membr,self.idx_membr) = self.membr.damping();
ele_mat(self.idx_plate,self.idx_plate) = self.plate.damping();
end
%% calculation of element force vector
function ele_vec = force(self)
%
ele_vec = zeros(20,1);
%
ele_vec(self.idx_membr) = self.membr.force();
ele_vec(self.idx_plate) = self.plate.force();
end
%------------------------------------------------------------------
%% post-process (compute internal forces)
function postprocess(self, u, du, ddu)
if ~exist('ddu', 'var')
ddu = zeros(length(u),1);
end
if ~exist('du', 'var')
du = zeros(length(u),1);
end
self.membr.postprocess(u(self.idx_membr), du(self.idx_membr), ddu(self.idx_membr));
self.plate.postprocess(u(self.idx_plate), du(self.idx_plate), ddu(self.idx_plate));
self.internal = [ self.membr.internal * self.sect.t ; self.plate.internal ];
end
%------------------------------------------------------------------
%% graphical output: plot element
function plot(self, config, scale)
% my nodes
n1 = self.node(1);
n2 = self.node(2);
n3 = self.node(3);
n4 = self.node(4);
% coordinates and displacements
x1 = [ n1.ref(1), n2.ref(1), n3.ref(1), n4.ref(1) ];
x2 = [ n1.ref(2), n2.ref(2), n3.ref(2), n4.ref(2) ];
w1 = [ n1.lhs( find(n1.dof == mafe.DofType.Disp1) ),
n2.lhs( find(n2.dof == mafe.DofType.Disp1) ),
n3.lhs( find(n3.dof == mafe.DofType.Disp1) ),
n4.lhs( find(n4.dof == mafe.DofType.Disp1) ) ]';
w2 = [ n1.lhs( find(n1.dof == mafe.DofType.Disp2) ),
n2.lhs( find(n2.dof == mafe.DofType.Disp2) ),
n3.lhs( find(n3.dof == mafe.DofType.Disp2) ),
n4.lhs( find(n4.dof == mafe.DofType.Disp2) ) ]';
w3 = [ n1.lhs( find(n1.dof == mafe.DofType.Disp3) ),
n2.lhs( find(n2.dof == mafe.DofType.Disp3) ),
n3.lhs( find(n3.dof == mafe.DofType.Disp3) ),
n4.lhs( find(n4.dof == mafe.DofType.Disp3) ) ]';
%
% defaults
lcolor = [0.6 0.6 0.6];
mcolor = [0.2 0.2 0.2];
lwidth = 0.5;
mwidth = 10.;
ecolor = [0.4 0.4 0.4];
marker = '.';
%
switch config
case 'deformed'
XX1 = x1 + w1*scale;
XX2 = x2 + w2*scale;
XX3 = w3*scale;
lcolor = [43.9, 50.6, 0.0]/100;
mcolor = [43.9, 50.6, 0.0]/100;
falpha = 0.2;
lwidth = 1.0;
mwidth = 20.;
case 'normal11'
XX1 = x1;
XX2 = x2;
lcolor = self.internal(1,:);
falpha = 1.0;
ecolor = 'none';
marker = 'none';
case 'normal22'
XX1 = x1;
XX2 = x2;
lcolor = self.internal(2,:);
falpha = 1.0;
ecolor = 'none';
marker = 'none';
case 'normal12'
XX1 = x1;
XX2 = x2;
lcolor = self.internal(3,:);
falpha = 1.0;
ecolor = 'none';
marker = 'none';
case 'moment11'
XX1 = x1;
XX2 = x2;
XX3 = w3*0.0;
lcolor = self.internal(4,:);
falpha = 1.0;
ecolor = 'none';
marker = 'none';
case 'moment22'
XX1 = x1;
XX2 = x2;
XX3 = w3*0.0;
lcolor = self.internal(5,:);
falpha = 1.0;
ecolor = 'none';
marker = 'none';
case 'moment12'
XX1 = x1;
XX2 = x2;
XX3 = w3*0.0;
lcolor = self.internal(6,:);
falpha = 1.0;
ecolor = 'none';
marker = 'none';
case 'vonmises'
XX1 = x1;
XX2 = x2;
XX3 = w3*0.0;
% calculate nodal von mises stress at upper surface (z=+t/2) from moment tensor
factor = 6. / self.sect.t^2;
s11 = self.internal(4,:) * factor;
s22 = self.internal(5,:) * factor;
s12 = self.internal(6,:) * factor;
sigma_v = sqrt( s11.^2 + s22.^2 - s11.*s22 + 3*s12.^2 );
lcolor = sigma_v;
falpha = 1.0;
ecolor = 'none';
marker = 'none';
case 'tensor'
XX1 = x1;
XX2 = x2;
XX3 = w3*0.0;
lcolor = [1.0 1.0 1.0];
mcolor = [0.1 0.1 0.1];
falpha = 1.0;
lwidth = 0.01;
mwidth = 0.01;
case 'shear13'
XX1 = x1;
XX2 = x2;
XX3 = w3*0.0;
lcolor = self.internal(7,:);
falpha = 1.0;
ecolor = 'none';
marker = 'none';
case 'shear23'
XX1 = x1;
XX2 = x2;
XX3 = w3*0.0;
lcolor = self.internal(8,:);
falpha = 1.0;
ecolor = 'none';
marker = 'none';
otherwise
XX1 = x1;
XX2 = x2;
XX3 = w3*0.0;
lcolor = [0.6 0.6 0.6];
mcolor = [0.2 0.2 0.2];
falpha = 1.0;
lwidth = 0.5;
mwidth = 12.;
end
patch(XX1, XX2, XX3, lcolor, ...
'EdgeColor', ecolor, ...
'LineWidth', lwidth, ...
'Marker', marker, ...
'MarkerSize', mwidth, ...
'FaceAlpha', falpha);
% put a tensor representation if required
if strcmp(config, 'tensor')
% compute tensor representation
S = [ mean(self.internal(4,:)), mean(self.internal(6,:));
mean(self.internal(6,:)), mean(self.internal(5,:)) ];
[D,M] = eig(S);
M = M * scale;
% compute center point of quadrilateral to place tensor cross
C = self.N([0,0]); CX1 = C*XX1'; CX2 = C*XX2';
% plot tensor cross
if M(1,1) < 0.0
color1 = 'b';
else
color1 = 'r';
end
if M(2,2) < 0.0
color2 = 'b';
else
color2 = 'r';
end
%plot( CX1, CX2, '*k' )
plot( CX1+M(1,1)*[-D(1,1) +D(1,1)], CX2+M(1,1)*[-D(2,1) +D(2,1)], color1, 'linewidth', 2.0 );
plot( CX1+M(2,2)*[-D(1,2) +D(1,2)], CX2+M(2,2)*[-D(2,2) +D(2,2)], color2, 'linewidth', 2.0 );
end
end
%% print element data
function str = print(self)
%
str_s11 = sprintf('%+4.3e ', self.internal(1,:)); % nodal s11 stresses
str_s22 = sprintf('%+4.3e ', self.internal(2,:)); % nodal s22 stresses
str_s12 = sprintf('%+4.3e ', self.internal(3,:)); % nodal s12 stresses
%
str_m11 = sprintf('%+4.3e ', self.internal(4,:)); % nodal s11 moments
str_m22 = sprintf('%+4.3e ', self.internal(5,:)); % nodal s22 moments
str_m12 = sprintf('%+4.3e ', self.internal(6,:)); % nodal s12 moments
%
str_q13 = sprintf('%+4.3e ', self.internal(7,:)); % nodal q13 shear forces
str_q23 = sprintf('%+4.3e ', self.internal(8,:)); % nodal q23 shear forces
str_1 = [ '\n normal: ' str_s11 '| ' str_s22 '| ' str_s12 ];
str_2 = [ '\n shear : ' str_q13 '| ' str_q23 ];
str_3 = [ '\n moment: ' str_m11 '| ' str_m22 '| ' str_m12 ];
str = [ str_1 str_2 str_3 ];
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