% __________________________________________________________________________
%| 64NODE-192DOF HEXAHEDRAL SOLID ISOPARAMETRIC F.E.M. [A.] (C)(R)|
%|########################################## |
%| MECHANISM | Corner upperside 4X(1Node-3Dof)[ux,vy,wz] |
%| | middle 8X(1Node-3Dof)[ux,vy,wz] |
%| | underside 4X(1Node-3Dof)[ux,vy,wz] |
%| | Side upperside 8X(1Node-3Dof)[ux,vy,wz] |
%| | middle 8X(1Node-3Dof)[ux,vy,wz] |
%| | underside 8X(1Node-3Dof)[ux,vy,wz] |
%| | Centroidal upperside 4X(1Node-3Dof)[ux,vy,wz] |
%| | middle 8X(1Node-3Dof)[ux,vy,wz] |
%| | underside 4X(1Node-3Dof)[ux,vy,wz] |
%|---------------------|----------------------------------------------------|
%| ANALYSIS | Shape function:Parametric (Pascal tree) |
%| | Stifness :Topology+Accumulate method |
%| | Solve Equation:Cholesky-[L][D][u] |
%|--------------------------------------------------------------------------|
%|SUB FUNCTIONS |
%|1- AconnectH64(pxn,pyn,pzn,e,n,J) [AconnectH64.m] |
%|2- threedimensionalstress(sigmax,sigmay,sigmaz,tauxy,tauyz,tauxz,1) |
%|--------------------------------------------------------------------------|
%| Maxima 5.9.0.9beta2 http://maxima.sourceforge.net |
%| Using Lisp Kyoto Common Lisp GCL 2.6.3 (aka GCL) |
%| Distributed under the GNU Public License. See the file COPYING. |
%| Dedicated to the memory of William Schelter. |
%| This is a development version of Maxima. The function bug_report() |
%| provides bug reporting information. |
%|__________________________________________________________________________|
%Shape function terms (Pascal tree)
% 1
% x,y,z
% x^2,y^2,z^2,x*y,y*z,x*z
% x^3,y^3,z^3,x^2*y,x*y^2,y*z^2,y^2*z,x*z^2,x^2*z,x*y*z
% x^3*y,x^2*y^2,x*y^3,y^3*z,y^2*z^2,y*z^3,x*z^3,x^2*z^2,x^3*z,x^2*y*z,x*y^2*z,x*y*z^2
% x^3*y^2,x^2*y^3,y^3*z^2,y^2*z^3,x^2*z^3,x^3*z^2,x^3*y*z,x*y^3*z,x*y*z^3,x^2*y*z^2,x*y^2*z^2
% x^3*y^3,y^3*z^3,x^3*z^3,x^3*y^2*z,x^2*y^3*z,x*y^3*z^2,x*y^2*z^3,x^2*y*z^3,x^3*y*z^2,x^2*y^2*z,x^2*y^2*z^2
% x^3*y^3*z,x*y^3*z^3,x^3*y*z^3,x^3*y^2*z^2,x^2*y^3*z^2,x^2*y^2*z^3
% x^3*y^3*z^2,x^2*y^3*z^3,x^3*y^2*z^3 %x^3*y^3*z^3
%ELEMENT PARAMETRIC SHAPE FUNCTIONS
%N1 : 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N2 : 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N3 : 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N4 : 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N49: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N50: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N51: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N52: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N5 : 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N6 : 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N9 : 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N10: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N53: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N54: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N57: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N58: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N7 : 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N8 : 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N11: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N12: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N55: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N56: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N59: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N60: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N17: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N18: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N19: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N20: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N33: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N34: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N35: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N36: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N21: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N22: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N25: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N26: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N37: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N38: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*1/16*(1-2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N41: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N42: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*1/16*(1+2*y/b)*(9*(2*y/b)^2-1)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N23: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N24: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N27: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N28: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N39: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N40: 1/16*(1+2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N43: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N44: 1/16*(1-2*x/a)*(9*(2*x/a)^2-1)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N29: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N30: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N31: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N32: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*9/16*(1-(2*z/c)^2)*(1-3*2*z/c);
%N45: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N46: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N47: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N48: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*9/16*(1-(2*z/c)^2)*(1+3*2*z/c);
%N13: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N14: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N15: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N16: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*1/16*(1-2*z/c)*(9*(2*z/c)^2-1);
%N61: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N62: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*9/16*(1-(2*y/b)^2)*(1-3*2*y/b)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N63: 9/16*(1-(2*x/a)^2)*(1+3*2*x/a)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
%N64: 9/16*(1-(2*x/a)^2)*(1-3*2*x/a)*9/16*(1-(2*y/b)^2)*(1+3*2*y/b)*1/16*(1+2*z/c)*(9*(2*z/c)^2-1);
clc
clear
%##################INPUT DATA###############
%####################### MATERIALS PROPERTIES
E=2E6; %shell material elasticity constant;
v=0.3; %shell material poission ratio;
%#######################
%####################### AUTOMATIC MESH DATA
totala = 2.00; %"x" axis direction total length (m)
totalb = 4.00; %"y" axis direction total width (m)
totalc = 1.00; %"z" axis direction total heigth (m)
meshx = 2; %"x" axis direction mesh number
meshy = 2; %"y" axis direction mesh number
meshz = 1; %"z" axis direction segment value
segment = 2; %Stress "z" axis directions values
%#######################
gridspacex = totala/meshx; % along to "x" axis direction mesh space
gridspacey = totalb/meshy; % along to "y" axis direction mesh space
gridspacez = totalc/meshz; % along to "z" axis direction mesh space
%Solid element constant dimensions
a=gridspacex ;
b=gridspacey ;
c=gridspacez ;
%Mesh error control
if meshx <1 || meshy <1 || meshz <1;% <>-----------(<>)
display('Open dimension system mesh and this mesh value > 1')% |
% [meshx meshy] |
error('mesh values')% |
else if meshx >10 || meshy>10 || meshz >10; % [>-----------(<>)
display('system analysis mesh is big')% |
end% |
end% <>-------------|
%_____________________________________GLOBAL COORDINATES
dim=1;
for ars=1:3*meshz+1 %>>-------------------------|
for sut=1:meshy*3+1 %>>-------------------| | |
for sat=1:meshx*3+1 %>>----------------| | | |
% [mod(sat,3) mod(sut,3)]; % | | | |
x=(sat-1)*gridspacex/3; % | | | |
y=(sut-1)*gridspacey/3; % | | | |
z=(ars-1)*gridspacez/3-totalc/2; % | | | |
Cor(dim,:)=[x y z]; % | | | |
dim=dim+1; % | | | |
end %<<-------sat------| | | |
end %<<-------sut---------| | |
end %<<-------------------------|
Number=size(Cor);
Node=Number(1); %Elements Number
%meshx=1;
%meshy=1;
%____________________________________POSITION MATRIX
no=1; %Pos(Element No,:)=[Element Node Number]
for k=1:3*meshz+1 % >>-------------------------|
for i=1:meshy;% >>------------------| | |
for j=1:meshx;% >>--------------| | | |
%Corner nodes
Ps(no,1)= (3*(3*meshx+1))*(i-1)+3*j-2 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,2)= (3*(3*meshx+1))*(i-1)+3*j+1 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,3)= (3*(3*meshx+1))*(i)+3*j+1 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,4)= (3*(3*meshx+1))*(i)+3*j-2 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
%Edge side nodes
%
%|__(5)__(6)___|
Ps(no,5)= (3*(3*meshx+1))*(i-1)+3*j-1 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,6)= (3*(3*meshx+1))*(i-1)+3*j ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,7)= (3*meshx+1)*(i)+2*(3*meshx+1)*(i-1)+3*j+1 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,8)= 2*(3*meshx+1)*(i)+(3*meshx+1)*(i-1)+3*j+1 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,9)= (3*(3*meshx+1))*(i)+3*j ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,10)= (3*(3*meshx+1))*(i)+3*j-1 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,11)= 2*(3*meshx+1)*(i)+(3*meshx+1)*(i-1)+3*j-2 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,12)= (3*meshx+1)*(i)+2*(3*meshx+1)*(i-1)+3*j-2 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
%Centroidal nodes
Ps(no,13)= (3*meshx+1)*(i)+2*(3*meshx+1)*(i-1)+3*j-1 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,14)= (3*meshx+1)*(i)+2*(3*meshx+1)*(i-1)+3*j ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,15)= 2*(3*meshx+1)*(i)+(3*meshx+1)*(i-1)+3*j ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
Ps(no,16)= 2*(3*meshx+1)*(i)+(3*meshx+1)*(i-1)+3*j-1 ...
+(k-1)*(3*meshx+1)*(3*meshy+1);
% | | |
% | | |
no=no+1 ;% | | |
end% <<-------j----------| | |
end% <<-------i-------------| |
end% <<------for(k)---------------|
for i=1:meshx*meshy%Total element number value
Pos(i,:)=[Ps(i,:) Ps(i+meshx*meshy,:) Ps(i+2*meshx*meshy,:) Ps(i+3*meshx*meshy,:)];
end
No=size(Pos,1);%Total element number
%Element`s Node Topology
for i=1:Node;
Re(i,:)=[1 1 1];
end
%Re(Node number,:)=[ux vy wz theta beta]
%############################SYSTEM SUPPORTS
Re(1,:) =[0 0 0]; %Simply supported
Re(7,:) =[0 0 0];
Re(43,:)=[0 0 0];
Re(49,:)=[0 0 0];
%#############################
Topology=Re;
%_________________________________ACCUMULATE
Nom=size(Re,2); %Element node d.o.f Number
sayman=0;
for i=1:Node;
for j=1:Nom;
if Re(i,j)==1 ;
sayman = sayman +1;
Re(i,j) = sayman;
end
end
end%_________________________________|
Item=sayman;
%_________________________________MODAL DISPLACEMENTS
for i=1:No
R(i,:)=[Re(Pos(i,1),:) Re(Pos(i,2),:) Re(Pos(i,3),:) Re(Pos(i,4),:) ...
Re(Pos(i,5),:) Re(Pos(i,6),:) Re(Pos(i,7),:) Re(Pos(i,8),:) ...
Re(Pos(i,9),:) Re(Pos(i,10),:) Re(Pos(i,11),:) Re(Pos(i,12),:) ...
Re(Pos(i,13),:) Re(Pos(i,14),:) Re(Pos(i,15),:) Re(Pos(i,16),:) ...
Re(Pos(i,17),:) Re(Pos(i,18),:) Re(Pos(i,19),:) Re(Pos(i,20),:) ...
Re(Pos(i,21),:) Re(Pos(i,22),:) Re(Pos(i,23),:) Re(Pos(i,24),:) ...
Re(Pos(i,25),:) Re(Pos(i,26),:) Re(Pos(i,27),:) Re(Pos(i,28),:) ...
Re(Pos(i,29),:) Re(Pos(i,30),:) Re(Pos(i,31),:) Re(Pos(i,32),:) ...
Re(Pos(i,33),:) Re(Pos(i,34),:) Re(Pos(i,35),:) Re(Pos(i,36),:) ...
Re(Pos(i,37),:) Re(Pos(i,38),:) Re(Pos(i,39),:) Re(Pos(i,40),:) ...
Re(Pos(i,41),:) Re(Pos(i,42),:) Re(Pos(i,43),:) Re(Pos(i,44),:) ...
Re(Pos(i,45),:) Re(Pos(i,46),:) Re(Pos(i,47),:) Re(Pos(i,48),:) ...
Re(Pos(i,49),:) Re(Pos(i,50),:) Re(Pos(i,51),:) Re(Pos(i,52),:) ...
Re(Pos(i,53),:) Re(Pos(i,54),:) Re(Pos(i,55),:) Re(Pos(i,56),:) ...
Re(Pos(i,57),:) Re(Pos(i,58),:) Re(Pos(i,59),:) Re(Pos(i,60),:) ...
Re(Pos(i,61),:) Re(Pos(i,62),:) Re(Pos(i,63),:) Re(Pos(i,64),:)];
end
%######################SYSTEM LOADS
%P(Re(Node number, Freedoom number))=[u v w theta beta]
P(Item)=0;
P(Re(4 +2*49,3))= -5;
P(Re(22+2*49,3))= -5;
P(Re(28+2*49,3))= -5;
P(Re(46+2*49,3))= -5;
%#######################
%"ELEMENT ELASTICITY MATRIX [C]";
C=E/((1+v)*(1-2*v))*[1-v , v , v , 0 , 0 , 0 ;
v ,1-v , v , 0 , 0 , 0 ;
v , v ,1-v, 0 , 0 , 0 ;
0 , 0 , 0 ,(1-2*v)/2, 0 , 0 ;
0 , 0 , 0 , 0 ,(1-2*v)/2, 0 ;
0 , 0 , 0 , 0 , 0 ,(1-2*v)/2];
%_____________________________GAUSS LEGENDRE CONSTANTS
%format long;
H=[0.171324492379170,0.360761573048139,0.467913934572691,0.467913934572691,0.360761573048139,0.171324492379170] ;
W=[0.932469514203152;
0.661209386466265;
0.238619186083197;
-0.238619186083197;
-0.661209386466265;
-0.932469514203152];
%_________________________________________________|
K(192,192,No)=0;
%_____________________________LOCAL AXIS ELEMENT STIFFNESS MATRIX
for s=1:No;%total element number
for u=1:64
pxn(u,:) =Cor(Pos(s,u),1); %x(u)
pyn(u,:) =Cor(Pos(s,u),2); %y(u)
pzn(u,:) =Cor(Pos(s,u),3); %z(u)
end
for k=1:6;
for j=1:6;
for i=1:6;
%[i j k s] Numerical integration numerator
%Parametric parameters is connection homogen parameters
%[i j k s] Numerical integration numerator
%Parametric parameters is connection homogen parameters
e=W(i);
n=W(j);
J=W(k);
%____FUNCTION____:Element shape functions partial derivatives.
[Jacobi,Hx,Hy,Hz]=AconnectH64(pxn,pyn,pzn,e,n,J);
%[A]=[Equbilibrium differantial matrix].[shell shape matrix]
%[A]=[d].[N]
%Element connection matrix
for f=1:64;
A(:,3*(f-1)+1,s)= [Hx(f); 0 ; 0 ; Hy(f) ; 0 ; Hz(f)];
A(:,3*(f-1)+2,s)= [ 0 ; Hy(f) ; 0 ; Hx(f) ; Hz(f); 0 ];
A(:,3*(f-1)+3,s)= [ 0 ; 0 ; Hz(f); 0 ; Hy(f); Hx(f)];
end
K(:,:,s)=K(:,:,s)+det(Jacobi)*H(i)*H(j)*H(k)*A(:,:,s)'*C*A(:,:,s);
end
end
end
fprintf ('Element integrate number = % .2f \n',s);
end%______________________________________________________________|
Ksis(Item,Item)=0;
%___________________________SYSTEM STIFFNESS MATRIX
for n=1:No;
for sat=1:192;
for sut=1:192;
if (R(n,sat)~=0)
if (R(n,sut)~=0);
Ksis(R(n,sut),R(n,sat))=Ksis(R(n,sut),R(n,sat)) + K(sat,sut,n);
end%
end%
end%
end%
end%_______________________________________________|
%Ksis:Global system stiffness matrix
clear n sat sut
%Ksis:System global stiffness matrix
%______________________________SYSTEM DISPLACEMENT
equation=size(Ksis);
if equation(1)~=rank(Ksis)
display('This system stiffness matrix is badly scaled')
% R
error('Control system support boundary conditions')
else
Ku = inv(Ksis);
D = Ku *P';
end%__________________________________________________|
clear equation Ku
%_________________________________MODAL DISPLACEMENT
for s = 1 : No;
for m = 1 :192;
u = R(s, m);
if u ~=0; Hu(s,m,:) = D(u) ;
else Hu(s,m,:) = 0 ;
end%
end%
end%__________________________________________________|
%Hu:per element node`s displacements
clear s m u
%Information
%This stress connection the global system axis components (x,y,z) and
%element stress infinite values depends this (x,y,z). This region
%per element has edge-middle node position stress values calculated.
fig=1;%Figure value
%__________________________________________________STRESS AND STRAIN ANALYSIS
for J=-c/2:c/segment:+c/2 ;%"z" axis direction position
for s=1:No;
for i=1:16;%Element stress value node numbers
switch(i)
case {1},
% x=-a/2; y=-b/2; %Element node for (1)
e=-1 ; n=-1 ; %Homogen coordinates
case {2},
% x=a/2 ; y=-b/2; %Element node for (2)
e= 1 ; n=-1 ; %Homogen coordinates
case {3},
% x=a/2;y=b/2; %Element node for (3)
e= 1 ; n= 1 ; %Homogen coordinates
case {4},
% x=-a/2;y=b/2; %Element node for (4)
e=-1 ; n=1 ; %Homogen coordinates
case {5},
e=-1/3 ; n=-1 ; %Homogen coordinates(5)
case {6},
e= 1/3 ; n=-1 ; %Homogen coordinates(6)
case {7},
e= 1 ; n=-1/3 ;%Homogen coordinates(7)
case {8},
e= 1 ; n= 1/3 ;%Homogen coordinates(8)
case {9},
e= 1/3 ; n= 1 ; %Homogen coordinates
case {10},
e= -1/3; n= 1 ; %Homogen coordinates
case {11},
e= -1 ; n=1/3 ; %Homogen coordinates
case {12},
e= -1 ; n=-1/3 ;%Homogen coordinates
case {13},
e=-1/3 ; n=-1/3 ;%Homogen coordinates
case {14},
e= 1/3 ; n=-1/3 ;%Homogen coordinates
case {15},
e= 1/3 ; n=1/3 ; %Homogen coordinates
case {16},
e=-1/3 ; n=1/3 ; %Homogen coordinates(16)
end
% x=gridspacex/2*e;%parametric axis convert homogen axis
% y=gridspacey/2*n;
%____FUNCTION____:Element shape functions partial derivatives.
[Jacobi,Hx,Hy,Hz]=AconnectH64(pxn,pyn,pzn,e,n,J);
%Element Connection Matrix [A]
for f=1:64;
A(:,3*(f-1)+1,s)= [Hx(f); 0 ; 0 ; Hy(f) ; 0 ; Hz(f)];
A(:,3*(f-1)+2,s)= [ 0 ; Hy(f) ; 0 ; Hx(f) ; Hz(f); 0 ];
A(:,3*(f-1)+3,s)= [ 0 ; 0 ; Hz(f); 0 ; Hy(f); Hx(f)];
end
Qeglobal(:,i,s)=A(:,:,s)*Hu(s,:)'; %Element global node`s Strains
Qsglobal(:,i,s)=C*Qeglobal(:,i,s); %Element global node`s Stresses
end
end
fprintf (' ###### SOLID "z" axis direction SEGMENT##### z=[% 2.4f]m \n',z);
%____________________________GLOBAL STRAINS
display(' PER ELEMENTS NODE GLOBAL STRAINS ===[ex ey ez xy yz xz]global');
for s=1:No;
fprintf(' Element number = ##%.f## \n',s);
fprintf(' Node number [1] [2] [3] [4] \n');
fprintf(' [% .7f] [% .7f] [% .7f] [% .7f] \n',Qeglobal(:,:,s));
fprintf(' \n');
end%______________________________|
fprintf(' \n \n \n')
%Qeglobal
%____________________________GLOBAL STRESSES
display (' PER ELEMENTS NODE GLOBAL STRESSES===[ox oy oz Txy Tyz Txz]global');
for s=1:No;
fprintf(' Element number = ##%.f## \n',s);
fprintf(' Node number [1] [2] [3] [4] \n');
fprintf(' [% .7f] [% .7f] [% .7f] [% .7f] \n',Qsglobal(:,:,s));
fprintf(' \n');
end%_______________________________|
fprintf(' \n \n \n')
%############# GRAPHICAL INTERFACE
px=1:3*meshx+1;%meshgridx
py=1:3*meshy+1;%meshgridy
for n=1:No
for i=1:16 %Element underside nodes
V1(Pos(n,i),:)=Qsglobal(1,i,n);%sigmaQX
V2(Pos(n,i),:)=Qsglobal(2,i,n);%sigmaQy
V3(Pos(n,i),:)=Qsglobal(3,i,n);%sigmaQz
V4(Pos(n,i),:)=Qsglobal(4,i,n);%tau Txy
V5(Pos(n,i),:)=Qsglobal(5,i,n);%tau Tyz
V6(Pos(n,i),:)=Qsglobal(6,i,n);%tau Txz
end
end
%Buradan eleman gerilme deerleri mesh fonksiyonuna aktarlyor
nm=1;
for i=1:3*meshx+1 ;
for j=1:3*meshy+1 ;
Qs1(i,j,:)=V1(nm) ;
Qs2(i,j,:)=V2(nm) ;
Qs3(i,j,:)=V3(nm) ;
Qs4(i,j,:)=V4(nm) ;
Qs5(i,j,:)=V5(nm) ;
Qs6(i,j,:)=V6(nm) ;
nm=nm+1;
end
end
tit=['Element stresses values of homogen axis direction (z)= ' num2str(J) ' metre'];
%Select graphical interface color map
%colormap(hsv)
%colormap(jet)
%colormap(cool)
%colormap(hot)
%colormap(gray)
%colormap(pink)
%colormap(copper)
figure(fig)
%figure(1)
subplot(2,3,1);contourf(Qs1);title('Sigmax Qx')
subplot(2,3,2);contourf(Qs2);title('Sigmay Qy')
subplot(2,3,3);contourf(Qs3);title('Sigmaz Qz')
subplot(2,3,4);contourf(Qs4);title('Tauxy Txy')
subplot(2,3,5);contourf(Qs5);title('Tauyz Tyz')
xlabel(tit);
subplot(2,3,6);contourf(Qs6);title('Tauxz Txz')
fig=fig+1;
%figure(fig)
%subplot(2,3,1);pcolor(px,py,Qs1);title('Sigmax Qx');shading interp;
%subplot(2,3,2);pcolor(px,py,Qs2);title('Sigmay Qy');shading interp;
%subplot(2,3,3);pcolor(px,py,Qs3);title('Sigmaz Qz');shading interp;
%subplot(2,3,4);pcolor(px,py,Qs4);title('Tauxy Txy');shading interp;
%subplot(2,3,5);pcolor(px,py,Qs5);title('Tauyz Tyz');shading interp;
%Xlabel(tit);
%subplot(2,3,6);pcolor(px,py,Qs6);title('Tauxz Qxz');shading interp;
%fig=fig+1;
end %_______________________________|
clear fig i j nm px py no tit x y z
%NODE REACTIONS
for s=1 :No;
Pg(:,s) = K(:,:,s)*Hu(s,:)'; %Global reactions
end
%_______________________________GLOBAL NODE DISPLACEMENT
display ('====PER ELEMENTS NODE GLOBAL DISPLACEMENT===[Du Dv Dw]global');
for s=1:size(Hu',1);
fprintf(' [% .7f] [% .7f] [% .7f] [% .7f] \n',Hu(:,s)');
%fprintf(' [% .7f]%.f [% .7f]%.f [% .7f]%.f [% .7f]%.f \n',Pg(s,1),R(1,s)',Pg(s,2),R(2,s)',Pg(s,3),R(3,s)',Pg(s,4),R(4,s)');
if mod(s,3)==0; fprintf(' \n');end;
end%____________________________________|
%_______________________________GLOBAL REACTIONS
display ('====PER ELEMENTS NODE GLOBAL REACTIONS===[Fx Fy Fz]global');
for s=1:size(Pg,1);
fprintf(' [% .7f] [% .7f] [% .7f] [% .7f] \n',Pg(s,:));
%fprintf(' [% .7f]%.f [% .7f]%.f [% .7f]%.f [% .7f]%.f \n',Pg(s,1),R(1,s)',Pg(s,2),R(2,s)',Pg(s,3),R(3,s)',Pg(s,4),R(4,s)');
if mod(s,3)==0; fprintf(' \n');end;
end%____________________________________|
clear s Node Nom
clear V1 V2 V3 V4 V5 V6 Qs1 Qs2 Qs3 Qs4 Qs5 Qs6
for i=1:No
for j=1:16
sigmax=Qsglobal(1,j,i);
sigmay=Qsglobal(2,j,i);
sigmaz=Qsglobal(3,j,i);
tauxy=Qsglobal(4,j,i);
tauyz=Qsglobal(5,j,i);
tauxz=Qsglobal(6,j,i);
[taumax,taumin,sigman1,sigman2,sigman3,sigmavonmises]= ...
threedimensionalstress(sigmax,sigmay,sigmaz,tauxy,tauyz,tauxz,1);
%[taumax,taumin,sigman1,sigman2,sigman3,sigmavonmises]= ...
%threedimensionalstress(sigmax,sigmay,sigmaz,tauxy,tauyz,tauxz,0)
Qss(:,j,i)=[taumax
taumin
sigman1
sigman2
sigman3
sigmavonmises];
end
end
clear No
%Ekstrem stress values
Qss
%################## INFORMATION
%Global axis direction per element stresses components
%Qsglobal(:,:,Element number)
%Qsglobal(:,:,100) =
%
% 1.0e+003 *
% [1] [2] [3] [4] Element node number--->
% 1.8487 1.4243 0.0098 0.4341 ox (Kn/m**2)
% 1.8487 0.4341 0.0098 1.4243 oy
% 0 0 0 0 oz
% -1.1700 -1.6651 -2.1602 -1.6651 Txy
% 0.0385 0 -0.0094 -0.0418 Tyz
% 0.0385 -0.0418 -0.0094 0 Txz
% /|\
% Stress components |
%
%Qss(:,:,4) =
% Node(1) Node(1) Node(1) Node(1)
% 6.2920 8.7557 7.5749 7.9609 taumax (Kn/m**2)
% -6.2920 -8.7557 -7.5749 -7.9609 taumin
% 10.4896 7.1814 -1.9850 3.8175 sigmamax1 (max)
% 0.4407 0.0157 -6.6046 -5.4689 sigmamax2 (med)
% -2.0944 -10.3300 -17.1349 -12.1044 sigmamax3 (min)
% 11.5275 15.2484 13.4489 13.8523 sigmamavonmises
% /|\
% Stress components |
