kakar133
Member level 2
- Joined
- Sep 16, 2010
- Messages
- 49
- Helped
- 0
- Reputation
- 0
- Reaction score
- 0
- Trophy points
- 1,286
- Location
- USA
- Activity points
- 1,579
hi every one I am trying calculate stiffness matrix for 2D fem and using quadrilateral element(2 elements). but i am facing problem with the integration.
it is giving me warnings with out calculating stiffness matrix. the warning is
Warning: Cannot attach the property of being close to the limit point to limit
variable. [limit]
please help me
syms z n;
epsilon=4*8.86e-12;
u=4*pi*1e-7;
w=2*pi*300e6;
si(1)=1/4*(1-z).*(1-n);
si(2)=1/4*(1+z).*(1-n);
si(3)=1/4*(1+z).*(1+n);
si(4)=1/4*(1-z).*(1+n);
for e=1:1:2;
J11(e) =1/4*(-1*(1-n)*x(e)+(1-n)*x(e+1)+(1+n)*x(e+2)-(1+n)*x(e+3));
J12(e) =1/4*(-1*(1-n)*y(e)+(1-n)*y(e+1)+(1+n)*y(e+2)-(1+n)*y(e+3));
J21(e) =1/4*(-1*(1-z)*x(e)-(1-z)*x(e+1)+(1+z)*x(e+2)+(1+z)*x(e+3));
J22(e) =1/4*(-1*(1-z)*y(e)-(1-z)*y(e+1)+(1+z)*y(e+2)+(1+z)*y(e+3));
j_d(e)=J11(e)*J22(e)-J12(e)*J21(e);
si_x(e,1)=(1/(4.*j_d(e))).*(-1*J22(e).*(1-n)+J12(e).*(1-z));
si_y(e,1)=(1/(4.*j_d(e))).*(J21(e).*(1-n)-J11(e).*(1-z));
si_x(e,2)=(1/(4.*j_d(e))).*(J22(e).*(1-n)+J12(e).*(1+z));
si_y(e,2)=(1/(4.*j_d(e))).*(-1*J21(e).*(1-n)-J11(e).*(1+z));
si_x(e,3)=(1/(4.*j_d(e))).*(J22(e).*(1+n)-J12(e).*(1+z));
si_y(e,3)=(1/(4.*j_d(e))).*(-1*J21(e).*(1+n)+J11(e).*(1+z));
si_x(e,4)=(1/(4.*j_d(e))).*(-1*J22(e).*(1+n)-J12(e).*(1-z));
si_y(e,4)=(1/(4.*j_d(e))).*(J21(e).*(1+n)+J11(e).*(1-z));
for j=1:1:4;
for i=1:1:4
K(j,i)= double(int(int((((1/epsilon)*(si_x(e,j).*si_x(e,i)+si_y(e,j).*si_y(e,i))-((w^2)*u*si(j).*si(i))).*j_d(e)),z,-1,1),n,-1,1));
%K(j,i)=dblquad(f,-1,1,-1,1);
end
end
end
it is giving me warnings with out calculating stiffness matrix. the warning is
Warning: Cannot attach the property of being close to the limit point to limit
variable. [limit]
please help me
syms z n;
epsilon=4*8.86e-12;
u=4*pi*1e-7;
w=2*pi*300e6;
si(1)=1/4*(1-z).*(1-n);
si(2)=1/4*(1+z).*(1-n);
si(3)=1/4*(1+z).*(1+n);
si(4)=1/4*(1-z).*(1+n);
for e=1:1:2;
J11(e) =1/4*(-1*(1-n)*x(e)+(1-n)*x(e+1)+(1+n)*x(e+2)-(1+n)*x(e+3));
J12(e) =1/4*(-1*(1-n)*y(e)+(1-n)*y(e+1)+(1+n)*y(e+2)-(1+n)*y(e+3));
J21(e) =1/4*(-1*(1-z)*x(e)-(1-z)*x(e+1)+(1+z)*x(e+2)+(1+z)*x(e+3));
J22(e) =1/4*(-1*(1-z)*y(e)-(1-z)*y(e+1)+(1+z)*y(e+2)+(1+z)*y(e+3));
j_d(e)=J11(e)*J22(e)-J12(e)*J21(e);
si_x(e,1)=(1/(4.*j_d(e))).*(-1*J22(e).*(1-n)+J12(e).*(1-z));
si_y(e,1)=(1/(4.*j_d(e))).*(J21(e).*(1-n)-J11(e).*(1-z));
si_x(e,2)=(1/(4.*j_d(e))).*(J22(e).*(1-n)+J12(e).*(1+z));
si_y(e,2)=(1/(4.*j_d(e))).*(-1*J21(e).*(1-n)-J11(e).*(1+z));
si_x(e,3)=(1/(4.*j_d(e))).*(J22(e).*(1+n)-J12(e).*(1+z));
si_y(e,3)=(1/(4.*j_d(e))).*(-1*J21(e).*(1+n)+J11(e).*(1+z));
si_x(e,4)=(1/(4.*j_d(e))).*(-1*J22(e).*(1+n)-J12(e).*(1-z));
si_y(e,4)=(1/(4.*j_d(e))).*(J21(e).*(1+n)+J11(e).*(1-z));
for j=1:1:4;
for i=1:1:4
K(j,i)= double(int(int((((1/epsilon)*(si_x(e,j).*si_x(e,i)+si_y(e,j).*si_y(e,i))-((w^2)*u*si(j).*si(i))).*j_d(e)),z,-1,1),n,-1,1));
%K(j,i)=dblquad(f,-1,1,-1,1);
end
end
end