index in position exceeds array bounds

%--------------------------BMRP COMPOSITE--------------%

% ----------------------nomenclature-------------------%

% l is the length of the given plate

% h is the width of given plate

% t is the total thickness

% nx is the total node in x direction

% ny is total nodes in y direction

% input paramteres

l=input('enter the length of plate');

h=input('enter the breadth of plate');

Nx=input('enter the no of elements in x direction');

Ny=input('enter the no of elements in Y direction');

t=input('enter the total thickness of laminate');

alpha=input('enter the correction factor');

T=input('enter the temperature');

C=input('enter the concentration');

T0=input('enter the reference temperature');

C0=input('enter the reference concentration');

R=input('enter the radius of curvature');

nne=8; % no of nodes per element

nodof=5; % degree of freedom per node

eldof=40;

dim=2;

nx=(2*Nx)+1;

ny=(2*Ny)+1;

nel=Nx*Ny; % total number of elements

nnd=341; % total number of nodes

p=nnd*nodof; % global degree of freedom

n=input('enter the number of plies');

z=-t/2:t/n:t/2;

z=z';

theta=[0;90;0;90;0;90;0;90;90;0;90;0;90;0;90;0];

gx=[-0.577;0.577;0.577;-0.577];

gy=[-0.577;-0.577;0.577;0.577];

nogp=4;

% material properties of lamina

E1=15.94E9;

E2=15.94E09;

G12=3.5709;

G13=G12;

G23=G12;

mu12=0.41;

mu21=0.41;

rho=1690;

alpha1=-0.3e-06;

alpha2=28.1e-06;

beta1=0;

beta2=0.44;

% generating mesh

[nodept,elnn]=rectangularmesh(l,h,nx,ny);

nodept=nodept';

for i=1:nel

k=8;

for j=3:5

t=elnn(i,j);

elnn(i,j)=elnn(i,k+1);

elnn(i,k+1)=t;

k=k-1;

end

end

% generating rigidity matrix or constitutive matrix D of laminate

[D,S1]=rigidity(n,theta,z,E1,E2,G12,G13,G23,mu12,mu21,alpha);

% generating the non mechanical force and moment for laminate

[Fn]=hygroforce(z,T0,C0,T,C,n,alpha1,alpha2,beta1,beta2,theta,S1);

% element stifness matrix and global stiffness matrix

% element mass matrix and global mass matrix

index=zeros(eldof,1);

kk=zeros(p,p); % global stiffness matrix

M=zeros(p,p); % global mass matrix

Knl=zeros(p,p); % global non linear stiffness matrix

for i=1:nel

kel=zeros(40);

mel=zeros(40);

Pel=zeros(40,1);

knl=zeros(40);

for r=1:nne

node(r)=elnn(i,nne+1);

end

[index]=elementdof(node,nne,nodof);

[coord]=cordinates(nne,nodept,elnn,dim,i); % coordinates of nodes of corresponding element

for j=1:nogp

zeeta1=gx(j);

eeta1=gy(j);

[J,derv]=jacobian(coord,zeeta1,eeta1);

d=det(J);

Jac1=inv(J);

der=Jac1*derv;

[B]=shapefunction(der,R,zeeta1,eeta1); % generating B matrix

kel=kel+(B'*D*B*d);

[N]=shapefunctionmatrix(zeeta1,eeta1);

[P]=densitymatrix(rho,z,n);

mel=mel+(N'*P*N*d);

Pel=Pel+(B'*Fn*d);

end

% extract system dofs associated with element

[kk]=globalstiff(index,eldof,kk,kel);

[M]=globalmass(index,eldof,M,mel);

delta=inv(kel)*Pel;

% Element non linear stifness matrix and assembling

for k=1:nogp

zeeta1=gx(k);

eeta1=gy(k);

[J,derv]=jacobian(coord,zeeta1,eeta1);

d=det(J);

Jac1=inv(J);

der=Jac1*derv;

[B]=shapefunction(der,R,zeeta1,eeta1);

Es=B*delta;

F=(D*Es)-Fn;

[S]=hygro(t,F);

[G]=nonlinear(zeeta1,eeta1,R,der);

knl=knl+(G'*S*G*d);

end

[Knl]=globalnonlinear(knl,Knl,index,eldof);

end

Thank you for your concern.I am sorry since i am new here i didnt think of it.This is the main program and i am getting error while generating index matrix of corresponding element.The required inputs are:All the function files are attached here.Please do respond

length=0.21 Temperature =300

breadth=0.21 concentration= 0

no of elemets in x drection =10 reference temperature =300

no of elements in y direction =10 reference concentration= 0

total thickness= 0.006 radius of curvature= 0.8673

correction factor= 0.8333 no of plies= 16

[KSSV]

Give the required inputs and copy the code here......with the given image snippet, without mentioning the inputs how do you think you can get help?

[Rohith Prakash]

I have updated everything.sorry for earlier inconvenience.Looking for your support

Answer 1

[madhan ravi]

nne-1