Materials and Computational Mechanics: 2011

Friday, December 9, 2011

programming

Hi Guys!

I found an interesting app on google chrome webstore which contains lots and lots of step by step lessons about programming languages, specifically C and C++. Its fascinating to know that there are videos, pdf files, cpp files and examples available on the website too. So if you are interested you can visit the website at:
website
The website describes itself as "WiBit.net is a video tutorial web site offering cutting edge programming and computer tutorials."
It is good to know that they use free compiler GCC and notepad++ as the text editor.

Saturday, December 3, 2011

J2-Hardening plasticity

Below you may find the constitutive driver for the J2-Plasticity model with isotropic hardening, considering associated flow rule, which means that the flow potential is the same as the yield surface and therefore the normal to the yield surface is identical to the normal to the plastic potential surface.

function [sig,eps,state_new,E_tang]=J2hardening(oldsig,oldeps,Deps,state_old,modpar)
% function [sig,eps,state_new,E_tang]=J2hardening(oldsig,oldeps,Deps,state_old,modpar)
%----------- Description ---------
% J2 elasto plasticity with isotropic hardening constitutive driver
% ------------ INPUT -------------
% oldsig = stress vector from the last increment [9X1]
% oldeps = strain vector from the last increment [9X1]
% Deps = strain increment [9X1]
% state_old = structure array with the specified fields and values
% could contain eps_e, eps_p, kappa, mu, etc.
% modpar = contains material constants
% ----------- OUTPUT -------------
% sig = stress vector from the current increment [9X1]
% eps = strain vector from the current increment [9X1]
% state_new = updated structured values
% E_tang = Tangent stiffness matrix
% Written by
% Shah, 02 December 2011

E=modpar.E;
nu=modpar.nu;
H=modpar.H;
sy=modpar.sy;
oldeps_p=state_old.eps_p;
oldkap=state_old.kap;
K=E/3/(1-2*nu);
G=E/2/(1+nu);
v1=[1 1 1 0 0 0 0 0 0]';
v2=[0 0 0 1 1 1 1 1 1]';
I=diag(v1+0.5*v2);
Idev=I-(1/3)*(v1*v1');
Ee=2*G*Idev+K*(v1*v1');
I2=diag(v1+v2);
I2dev=I2-(1/3)*(v1*v1');
sigT=oldsig+Ee*Deps;
sigDT=I2dev*sigT;
sigET=sqrt(1.5)*norm(sigDT);
phiT=sigET-sy-oldkap;
sig_m1=(sigT(1)+sigT(2)+sigT(3))/3;
eps=oldeps+Deps;
if phiT<=0
sig=sigT;
mu=0;
kap=oldkap;
eps_p=oldeps_p;
E_tang=Ee;
else
mu=phiT/(3*G+H);
c=1-mu*3*G/sigET;
sigD=sigDT*c;
sig_m=(sigT(1)+sigT(2)+sigT(3))/3;
sig=sigD+sig_m.*v1;
nuT=(1.5)*(sigDT/sigET);
kap=oldkap+H*mu;
E_tang=Ee-2*G*((2*G/(3*G+H)).*((oldkap+sy)/sigET)*nuT*nuT'+mu*3*G*Idev/sigET);
eps_p=oldeps_p+mu;
end
state_new.eps_p=eps_p;
state_new.J2=0;
state_new.I1=0;
state_new.kap=kap;
state_new.mu=mu;

Friday, December 2, 2011

GID

There is another mesh generator GID which is quite a powerful and complete mesh tool. This software is not free, but you can download it for free and use it with some limitations such as number of elements or nodes. Here I post the download link, give it a try, you won't be disappointed.
http://www.gidhome.com/download
There are versions for Windows, Mac and Linux for 32bit and 64bit.

Friday, November 18, 2011

Transfer your mesh from GMSH to Matlab

I believe examples are on of the best ways to learn, so I start from an example.

1- build your geometry in GMSH (you can easily do that by studying the first tutorial of the software)

2- Mesh the geometry by going to the Mesh tab and clicking on the 2D. This will build a triangular non structured mesh

Homework 8 solution

Homework 8 solution (click to download from Scribd)

HW#8

HW8 problem description (click to download from Scribd)

Sunday, November 13, 2011

GMSH

Here is the link to a very useful free mesh generator and post processor software (GMSH) which is distributed under the terms of the GNU General Public License (GPL). http://geuz.org/gmsh/ There are versions available for Win, Mac and Linux.
You can find tutorials in the tutorial folder of the downloaded file.

Below is an example of unstructured mesh for a circular surface with a squared hole inside.

Friday, November 11, 2011

4 node quadrilateral element for non symmetrical rectangles

Below you may find the Matlab function for 4 node quadrilateral element which is capable of meshing non symmetric rectangles. I will try to do the same for triangles.

function [Edof,Ex,Ey]=quadmesh2(xcorner,ycorner,nx,ny)

% function [Edof,Ex,Ey]=quadmesh2(xcorner,ycorner,nx,ny)
%----------- Description ---------
% Generates quadrilateral 4 node elements
% ------------ INPUT -------------
% xcorner=[lower left, lower right, upper right, upper left]
% ycorner=[lower left, lower right, upper right, upper left]
% nx= No. elements in x direction
% ny= No. elements in y direction
% 4 3
% o----------o
% | |
% | |
% | |
% o-------- o
% 1 2
% ----------- OUTPUT -------------
% Ex= element x coordinates matrix [nrelem X 4]
% Ey= element y coordinates matrix [nrelem X 4]
% Edof= element DOF [nrelem X 9]
% Written by
% Shah, 11 November 2011

ld=abs(xcorner(2)-xcorner(1));
lu=abs(xcorner(3)-xcorner(4));
hr=abs(ycorner(3)-ycorner(2));
hl=abs(ycorner(4)-ycorner(1));
for i=1:ny
for j=1:nx
dl=(i-1)*(nx+1)+j;
ul=dl+(nx+1);
dr=dl+1;
ur=ul+1;
l1=(i-1)*(lu-ld)/ny+ld; l2=i*(lu-ld)/ny+ld;
h1=(j-1)*(hr-hl)/nx+hl; h2=j*(hr-hl)/nx+hl;
x_cd=(i-1)*(xcorner(4)-xcorner(1))/ny+xcorner(1);
x_cu=i*(xcorner(4)-xcorner(1))/ny+xcorner(1);
y_cl=(j-1)*(ycorner(2)-ycorner(1))/nx+ycorner(1);
y_cr=j*(ycorner(2)-ycorner(1))/nx+ycorner(1);
xy_dl=[(j-1)*(l1/nx)+x_cd (i-1)*(h1/ny)+y_cl];
xy_dr=[j*(l1/nx)+x_cd (i-1)*(h2/ny)+y_cr];
xy_ul=[(j-1)*(l2/nx)+x_cu i*(h1/ny)+y_cl];
xy_ur=[j*(l2/nx)+x_cu i*(h2/ny)+y_cr];
el=(i-1)*nx+j;
Edof(el,:)=[el,2*dl-1,2*dl,2*dr-1,2*dr,2*ur-1,2*ur,2*ul-1,2*ul];
Ex(el,:)=[xy_dl(1) xy_dr(1) xy_ur(1) xy_ul(1)];
Ey(el,:)=[xy_dl(2) xy_dr(2) xy_ur(2) xy_ul(2)];
end
end

figure(1);
hold on;
for i=1:length(Edof(:,1));
plot(Ex(i,[1:end 1]),Ey(i,[1:end 1]));
end
axis equal

3 node symmetric triangular elements

Below you may find a Matlab function to generate 3 node symmetric triangular elements as I promised.
Note that nx and ny are the number of rectangular elements in each direction.

function [Edof,Ex,Ey]=symtrimesh(xcorner,ycorner,nx,ny)

% function [Edof,Ex,Ey]=symtrimesh(xcorner,ycorner,nx,ny)
%----------- Description ---------
% Generates triangular 3 node symmetric elements
% ------------ INPUT -------------
% xcorner=[lower left, lower right, upper right, upper left]
% ycorner=[lower left, lower right, upper right, upper left]
% nx= No. of rectangulat elements in x direction
% ny= No. of rectangulat elements in y direction
% 7 8 9
% o---------o---------o
% | + | + |
% | + | + |
% | + | + |
% | + | + |
% 4o---------5---------o6
% | + | + |
% | + | + |
% | + | + |
% | + | + |
% o---------o---------o
% 1 2 3
% ----------- OUTPUT -------------
% Ex= element x coordinates matrix [nrelem X 3]
% Ey= element y coordinates matrix [nrelem X 3]
% Edof= element DOF [nrelem X 7]
% Written by
% Shah, 11 November 2011

l=xcorner(2)-xcorner(1);

h=ycorner(3)-ycorner(2);

for i=1:ny
for j=1:nx

dl=(i-1)*(2*nx+1)+2*j-1;
dm=dl+1;
dr=dl+2;
ul=dl+(2*nx+1);
um=ul+1;
ur=ul+2;
xy_dl=[(j-1)*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];
xy_dr=[j*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];
xy_dm=(xy_dl+xy_dr)./2;
xy_ul=[(j-1)*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
xy_ur=[j*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
xy_um=(xy_ul+xy_ur)./2;
el=4*nx*(i-1)+4*(j-1)+1;
if floor(i/2)-i/2~=0
Edof(el,:)=[el,2*dl-1,2*dl,2*um-1,2*um,2*ul-1,2*ul];
Edof(el+1,:)=[el+1,2*dl-1,2*dl,2*dm-1,2*dm,2*um-1,2*um];
Edof(el+2,:)=[el+2,2*dm-1,2*dm,2*dr-1,2*dr,2*um-1,2*um];
Edof(el+3,:)=[el+3,2*dr-1,2*dr,2*ur-1,2*ur,2*um-1,2*um];
Ex(el,:)=[xy_dl(1) xy_um(1) xy_ul(1)];
Ex(el+1,:)=[xy_dl(1) xy_dm(1) xy_um(1)];
Ex(el+2,:)=[xy_dm(1) xy_dr(1) xy_um(1)];
Ex(el+3,:)=[xy_dr(1) xy_ur(1) xy_um(1)];
Ey(el,:)=[xy_dl(2) xy_um(2) xy_ul(2)];
Ey(el+1,:)=[xy_dl(2) xy_dm(2) xy_um(2)];
Ey(el+2,:)=[xy_dm(2) xy_dr(2) xy_um(2)];
Ey(el+3,:)=[xy_dr(2) xy_ur(2) xy_um(2)];
else
Edof(el,:)=[el,2*dl-1,2*dl,2*dm-1,2*dm,2*ul-1,2*ul];
Edof(el+1,:)=[el+1,2*dm-1,2*dm,2*um-1,2*um,2*ul-1,2*ul];
Edof(el+2,:)=[el+2,2*dm-1,2*dm,2*ur-1,2*ur,2*um-1,2*um];
Edof(el+3,:)=[el+3,2*dm-1,2*dm,2*dr-1,2*dr,2*ur-1,2*ur];
Ex(el,:)=[xy_dl(1) xy_dm(1) xy_ul(1)];
Ex(el+1,:)=[xy_dm(1) xy_um(1) xy_ul(1)];
Ex(el+2,:)=[xy_dm(1) xy_ur(1) xy_um(1)];
Ex(el+3,:)=[xy_dm(1) xy_dr(1) xy_ur(1)];
Ey(el,:)=[xy_dl(2) xy_dm(2) xy_ul(2)];
Ey(el+1,:)=[xy_dm(2) xy_um(2) xy_ul(2)];
Ey(el+2,:)=[xy_dm(2) xy_ur(2) xy_um(2)];
Ey(el+3,:)=[xy_dm(2) xy_dr(2) xy_ur(2)];

end
end
end

figure(1);

hold on;
for i=1:length(Edof(:,1));

plot(Ex(i,[1:end 1]),Ey(i,[1:end 1]));

end
axis equal

3 node triangular element

Below you may find a Matlab function to generate 3 node triangular elements. The output will not be symmetric, so I will try to write one for symmetric triangular elements soon.
Note that nx and ny are the number of rectangular elements in each direction.

function [Edof,Ex,Ey]=trimesh(xcorner,ycorner,nx,ny)

% function [Edof,Ex,Ey]=trimesh(xcorner,ycorner,nx,ny)
%----------- Description ---------
% Generates triangular 3 node elements
% ------------ INPUT -------------
% xcorner=[lower left, lower right, upper right, upper left]
% ycorner=[lower left, lower right, upper right, upper left]
% nx= No. elements in x direction
% ny= No. elements in y direction
% 3 4
% o---------o
% | + |
% | + |
% | + |
% | + |
% o---------o
% 1 2
% ----------- OUTPUT -------------
% Ex= element x coordinates matrix [nrelem X 3]
% Ey= element y coordinates matrix [nrelem X 3]
% Edof= element DOF [nrelem X 7]
% Written by
% Shah, 11 November 2011

l=xcorner(2)-xcorner(1);

h=ycorner(3)-ycorner(2);

for i=1:ny
for j=1:nx

dl=(i-1)*(nx+1)+j;

ul=dl+(nx+1);

dr=dl+1;

ur=ul+1;

xy_dl=[(j-1)*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];

xy_dr=[j*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];

xy_ul=[(j-1)*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];

xy_ur=[j*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];

el=2*(i-1)*nx+2*j-1;

Edof(el,:)=[el,2*dl-1,2*dl,2*dr-1,2*dr,2*ul-1,2*ul];

Edof(el+1,:)=[el+1,2*dr-1,2*dr,2*ur-1,2*ur,2*ul-1,2*ul];

Ex(el,:)=[xy_dl(1) xy_dr(1) xy_ul(1)];

Ex(el+1,:)=[xy_dr(1) xy_ur(1) xy_ul(1)];

Ey(el,:)=[xy_dl(2) xy_dr(2) xy_ul(2)];

Ey(el+1,:)=[xy_dr(2) xy_ur(2) xy_ul(2)];

end
end

figure(1);

hold on;
for i=1:length(Edof(:,1));

plot(Ex(i,[1:end 1]),Ey(i,[1:end 1]));

end
axis equal

Thursday, November 10, 2011

4 node quadrilateral element

Below you may find a Matlab function written by me to generate 4 node quadrilateral elements. It is very easy and straight forward.

function [Edof,Ex,Ey]=quadmesh(xcorner,ycorner,nx,ny)

% function [Edof,Ex,Ey]=quadmesh(xcorner,ycorner,nx,ny)
%----------- Description ---------
% Generates quadrilateral 4 node elements
% ------------ INPUT -------------
% xcorner=[lower left, lower right, upper right, upper left]
% ycorner=[lower left, lower right, upper right, upper left]
% nx= No. elements in x direction
% ny= No. elements in y direction
% 4 3
% o---------o
% | |
% | |
% | |
% o---------o
% 1 2
% ----------- OUTPUT -------------
% Ex= element x coordinates matrix [nrelem X 4]
% Ey= element y coordinates matrix [nrelem X 4]
% Edof= element DOF [nrelem X 9]
% Written by
% Shah, 10 November 2011

l=xcorner(2)-xcorner(1);
h=ycorner(3)-ycorner(2);
for i=1:ny
for j=1:nx
dl=(i-1)*(nx+1)+j;
ul=dl+(nx+1);
dr=dl+1;
ur=ul+1;
xy_dl=[(j-1)*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];
xy_dr=[j*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];
xy_ul=[(j-1)*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
xy_ur=[j*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
el=(i-1)*nx+j;
Edof(el,:)=[el,2*dl-1,2*dl,2*dr-1,2*dr,2*ur-1,2*ur,2*ul-1,2*ul];
Ex(el,:)=[xy_dl(1) xy_dr(1) xy_ur(1) xy_ul(1)];
Ey(el,:)=[xy_dl(2) xy_dr(2) xy_ur(2) xy_ul(2)];
end
end

figure(1);
hold on;
for i=1:length(Edof(:,1));
plot(Ex(i,[1:end 1]),Ey(i,[1:end 1]));
end
axis equal

Wednesday, November 9, 2011

Homework 7 solution

Homework 7 solution (click to download from Scribd)

HW#7

HW7 problem description (click to download from Scribd)

Friday, October 28, 2011

start of the project

From know on, I am going to post finite element codes on this blog. I will start from easy examples and will go deeper and deeper to reach very complex codes. I am using Matlab software at the moment but as the codes get more complex I might use Fortran instead.

Wednesday, April 27, 2011

RVE

In a little bit from know I am going to write about Representative Volume Element (RVE), and the different techniques used for homogenizing this heterogeneous element; such as Voigt (Taylor), Reuss(Hill) and etc. assumptions for subscale.
Some of the familiar names in this topic are: Hashin and Shtrikman, Mori and Tanka, Zohdi and Weiggers.

Wednesday, April 20, 2011

Hello World

I am going to write about constitutive modeling of materials, strong format and weak form of the governing laws, discretization and finite element modeling. My plan is to cover a wide range of material behavior such as elastic, elasto-plastic, viscoplastic, damage, porosity, hyperelasticity, finite strain formulation and etc.