No BSD License  

Highlights from
Two step Taylor Galerkin algorithm

image thumbnail

Two step Taylor Galerkin algorithm

by

 

18 Jun 2007 (Updated )

The two step Taylor Galerkin algorithm applied to nonlinear structural dynamics

...
function [Error_RF,DE] = ...
          Mesh_Error(NOC,NE,UG,B0,Young,Poisson)
%
% SPECIFICATION: This function serves to calculate some parameter to
%                control the errors 
%
% MAIN VARIABLES
%
% Error_RF:  Average rotational of F (deformation gradient) per element
%            It must be zero becase the rotational of a gradient is 0.
% NOC:       Nodal Coordinates
% NE:        Number of elements
% UG:        Unknowns in global coordinates      
% B:         Bij is the partial of Ni respect to Xj (Voigt matrix)
% Young,Poisson: elastic parameters

% Initialize
Error_RF=0;
DE=0;

for N=1:NE
   
   Node1 = NOC(N,1);
   Node2 = NOC(N,2);
   Node3 = NOC(N,3);
   Node4 = NOC(N,4);
   
   % Voigt matrix
   BE = B0(1:4,1:2,N);
   
   % Local unknowns
   [UE] = Local(UG,Node1,Node2,Node3,Node4);
   
   % Deformation gradient F

   F1(1,1)=UE(3);   F1(1,2)=UE(4);  % F Node1
   F1(2,1)=UE(5);   F1(2,2)=UE(6);   
   F2(1,1)=UE(9);   F2(1,2)=UE(10); % F Node2
   F2(2,1)=UE(11);  F2(2,2)=UE(12); 
   F3(1,1)=UE(15);  F3(1,2)=UE(16); % F Node3
   F3(2,1)=UE(17);  F3(2,2)=UE(18);   
   F4(1,1)=UE(21);  F4(1,2)=UE(22); % F Node4
   F4(2,1)=UE(23);  F4(2,2)=UE(24);
   
   F_Ave=(F1+F2+F3+F4)/4;
   
   % 1st Piola-Kirchhoff stress tensor
   Piola1 = Piola(Young,Poisson,F1);
   Piola2 = Piola(Young,Poisson,F2);
   Piola3 = Piola(Young,Poisson,F3);
   Piola4 = Piola(Young,Poisson,F4);
   Piola_Ave=(Piola1+Piola2+Piola3+Piola4)/4;
   
   % Divergence of 1st P-K
   [DP] = Div_P(Piola1,Piola2,Piola3,Piola4,BE);
   
   % Linear momentum   
   p1 = UE(1:2,1);
   p2 = UE(7:8,1);
   p3 = UE(13:14,1);
   p4 = UE(19:20,1);
   p_Ave=(p1+p2+p3+p4)/4;
   
   % Linear momentum gradient
   [Grad_p]=Grad_Lin_Mom(p1,p2,p3,p4,BE);      
   
   % Rotational error averaged
  [ErrorE]   = Error_Rot_E(F1,F2,F3,F4,BE);
   Error_RF  = Error_RF+(ErrorE/NE);
   
   % Total derivative of entropy 
  [EntropyE] = Entropy_E(Grad_p,DP,p_Ave,Piola_Ave);
   DE        = DE + EntropyE;    
   
end; % End of for N=1:NE

% End of function Error
   

% ----------- function Error_RotF_E --------------

function [ErrorE]=Error_Rot_E(F1,F2,F3,F4,B);
%
% SPECIFICATION: This Function serves to calculate the error 
%                of the vector of Unknowns 
% DESCRIPTION OF MAIN VARIABLES:
% 
% INPUT
%
% Fi:        Deformation gradient of node i
%             | Fxx  Fxy |
%             | Fyx  Fyy |
% dt:        Increment of ime step
% BE:         Bij is the partial of Ni respect to Xj (Voigt matrix)
%
% OUTPUT
%
% RotF:     Rotational of the deformation gradient F
%              | Fxy,x - Fxx,y |
%              | Fyy,x - Fyx,y |
      
  RotF=zeros(2,1);

  DxFxy = F1(1,2)*B(1,1) + F2(1,2)*B(2,1) +...
          F3(1,2)*B(3,1) + F4(1,2)*B(4,1);
  DyFxx = F1(1,1)*B(1,2) + F2(1,1)*B(2,2) +...
          F3(1,1)*B(3,2) + F4(1,1)*B(4,2);
  DxFyy = F1(2,2)*B(1,1) + F2(2,2)*B(2,1) +...
          F3(2,2)*B(3,1) + F4(2,2)*B(4,1);
  DyFyx = F1(2,1)*B(1,2) + F2(2,1)*B(2,2) +...
          F3(2,1)*B(3,2) + F4(2,1)*B(4,2);

  RotF(1,1)=DxFxy-DyFxx;
  RotF(2,1)=DxFyy-DyFyx;
   
  ErrorE=sqrt(RotF(1)*RotF(1)+RotF(2)*RotF(2));
   
% End of function Error_Rot_F_E

%------------------------  function Local  ---------------------------
%
function [UE]=Local(UG,Node1,Node2,Node3,Node4);
%
% SPECIFICATION: This function serves to obtain 
%                the local unknowns from global coordinates
%
% UG: Global vector
% UE: Local vector for the element


   UE=zeros(24,1);
   
   % Node 1
   Start1 = 6*(Node1-1)+1;
   End1   = Start1+5;
   UE(1:6,1)   = UG(Start1:End1,1);
   
   % Node 2
   Start2 = 6*(Node2-1)+1;
   End2   = Start2+5;
   UE(7:12,1)  = UG(Start2:End2,1);
  
   % Node 3
   Start3 = 6*(Node3-1)+1;
   End3   = Start3+5;
   UE(13:18,1) = UG(Start3:End3,1);
   
   % Node 4
   Start4 = 6*(Node4-1)+1;
   End4   = Start4+5;
   UE(19:24,1) = UG(Start4:End4,1);
    
% End of function Local


%------- function Grad_Lin_Mom -----------------
function [Grad_p]=Grad_Lin_Mom(p1,p2,p3,p4,B);
%
% SPECIFICATION: This Fluxction serves to calculate the gradient 
%                of the vector of Unknowns 
%
% DESCRIPTION OF MAIN VARIABLES:
% 
% INPUT
%
% pi=   linear momentum of node i
%            p=[ px
%                py ]
% dt:   Step of time
% B:    Bij is the partial of Ni respect to Xj (Voigt matrix)
%
% OUTPUT
%
% Grad_p:   Gradient of Unknowns
%            | px,x    px,y |
%            | py,x    py,y |

Grad_p=zeros(2,2);

% First column:  derivatives with respect to X
Grad_p(1:2,1)=p1(1:2)*B(1,1)+p2(1:2)*B(2,1)+p3(1:2)*B(3,1)+p4(1:2)*B(4,1);

% Second column: derivatives with respect to Y
Grad_p(1:2,2)=p1(1:2)*B(1,2)+p2(1:2)*B(2,2)+p3(1:2)*B(3,2)+p4(1:2)*B(4,2);
   
% End of Grad_Lin_Mom function


%------- function Div_P ----------------- 
%
function [DP] = Div_P(P1,P2,P3,P4,B)
%
% SPECIFICATION: This function serves to calculate the divergence 
%                of the 1st Piola Kirchhoff stress tensor 

DP=zeros(2,1);

DP=P1(1:2,1)*B(1,1)+P2(1:2,1)*B(2,1)+P3(1:2,1)*B(3,1)+P4(1:2,1)*B(4,1)+...
   P1(1:2,2)*B(1,2)+P2(1:2,2)*B(2,2)+P3(1:2,2)*B(3,2)+P4(1:2,2)*B(4,2);

% End of Div_P function

%------- function Entropy_E -----------------  
%
function [EntropyE] = Entropy_E(Grad_p,DP,p_Ave,Piola_Ave)
%
% SPECIFICATION: This function serves to calculate the variation of the 
%                entropy with respect to time dE/dt in each element

EntropyE=trace(transpose(Grad_p)*Piola_Ave)+transpose(p_Ave)*DP;

% End of Entropy_E function

Contact us