CALFEM / calfem-matlab

CALFEM - a finite element toolbox for MATLAB
MIT License
86 stars 36 forks source link

Having issues getting the joint ball boundary condition to work #6

Closed layzeees closed 2 years ago

layzeees commented 2 years ago

Here is my code along with the frame structure with Edof....

received-384836533421647

I am having a value for the moment (supposed to be null) at the join ball located at the top center of the structure (x=4.5m ; y=10m)

moment

i have my boudary conditions set as bc=[1 0; 2 0; 3 0; 4 0; 5 0; 6 0; 24 0];

i've tried bc=[1 0; 2 0; 3 0; 4 0; 5 0; 6 0];
and i got the exact same results which means ... me setting 24(ddl)=0 isn't doing anything.

%----------------------------------------------------------------
 echo on

%----- Topology -------------------------------------------------

 Edof=[1  1  2  3 7  8  9;
       2  4  5  6 16  17  18;
       3  7  8  9 10  11  12;
       4  10  11  12 13  14  15;
       5  13  14  15 16  17  18;
       6  7  8  9 19  20  21;
       7  16  17  18 25  26  27;
       8  19  20  21 22  23  24;
       9  22  23  24 25  26  27]      

%----- Stiffness matrix K and load vector f ---------------------

 K=zeros(27);   f=zeros(27,1);
 f(11)=-20e3;
 f(14)=-20e3;

%----- Element stiffness and element load matrices  -------------

 Eh=36e6;
 Ah=0.16;
 Ih=2.133e-3;

 Ev=70e6;
 Av=0.1;
 Iv=1.333e-3;

 eph=[Eh Ah Ih];
 epv=[Ev Av Iv];

 ex1=[0 0];       ex2=[9 9];         ex3=[0 3];   
 ey1=[0 5];      ey2=[0 5];      ey3=[5 5];     

 ex4=[3 6];      ex5=[6 9];       ex6=[0 0];
 ey4=[5 5];      ey5=[5 5];      ey6=[5 10];

 ex7=[9 9];       ex8=[0 4.5];    ex9=[4.5 9];    
 ey7=[5 10];       ey8=[10 10];    ey9=[10 10];

 eq1=[0 0];
 eq2=[0 0];
 eq3=[0 0];
 eq4=[0 0];
 eq5=[0 0];
 eq6=[0 0];
 eq7=[0 0];
 eq8=[0 -20e3];
 eq9=[0 -20e3];

 Ke1=beam2e(ex1,ey1,epv);
 Ke2=beam2e(ex2,ey2,epv);
 Ke3=beam2e(ex3,ey3,eph);
 Ke4=beam2e(ex4,ey4,eph);
 Ke5=beam2e(ex5,ey5,eph);
 Ke6=beam2e(ex6,ey6,epv);
 Ke7=beam2e(ex7,ey7,epv);

 [Ke8,fe8]=beam2e(ex8,ey8,eph,eq8);
 [Ke9,fe9]=beam2e(ex9,ey9,eph,eq9);

%----- Assemble Ke into K ---------------------------------------

 K=assem(Edof(1,:),K,Ke1);
 K=assem(Edof(2,:),K,Ke2);
 K=assem(Edof(3,:),K,Ke3);
 K=assem(Edof(4,:),K,Ke4);
 K=assem(Edof(5,:),K,Ke5);
 K=assem(Edof(6,:),K,Ke6);
 K=assem(Edof(7,:),K,Ke7);
 [K,f]=assem(Edof(8,:),K,Ke8,f,fe8);
 [K,f]=assem(Edof(9,:),K,Ke9,f,fe9);

%----- Solve the system of equations and compute reactions ------

 bc=[1 0; 2 0; 3 0; 4 0; 5 0; 24 0; 6 0];   
 [a,r]=solveq(K,f,bc)

%----- Section forces -------------------------------------------

 Ed=extract(Edof,a);

 es1=beam2s(ex1,ey1,epv,Ed(1,:),eq1,21) 
 es2=beam2s(ex2,ey2,epv,Ed(2,:),eq2,21) 
 es3=beam2s(ex3,ey3,eph,Ed(3,:),eq3,21)
 es4=beam2s(ex4,ey4,eph,Ed(4,:),eq4,21)
 es5=beam2s(ex5,ey5,eph,Ed(5,:),eq5,21)
 es6=beam2s(ex6,ey6,epv,Ed(6,:),eq6,21)
 es7=beam2s(ex7,ey7,epv,Ed(7,:),eq7,21)
 es8=beam2s(ex8,ey8,eph,Ed(8,:),eq8,21)
 es9=beam2s(ex9,ey9,eph,Ed(9,:),eq9,21)

 %----- Draw deformed frame ---------------------------------------

 figure(1)
 plotpar=[2 1 0];
 eldraw2(ex1,ey1,plotpar);
 eldraw2(ex2,ey2,plotpar);
 eldraw2(ex3,ey3,plotpar);
 eldraw2(ex4,ey4,plotpar);
 eldraw2(ex5,ey5,plotpar);
 eldraw2(ex6,ey6,plotpar);
 eldraw2(ex7,ey7,plotpar);
 eldraw2(ex8,ey8,plotpar);
 eldraw2(ex9,ey9,plotpar);
 sfac=scalfact2(ex7,ey7,Ed(7,:),0.1);
 plotpar=[1 2 1];
 eldisp2(ex1,ey1,Ed(1,:),plotpar,sfac);
 eldisp2(ex2,ey2,Ed(2,:),plotpar,sfac);
 eldisp2(ex3,ey3,Ed(3,:),plotpar,sfac);
 eldisp2(ex4,ey4,Ed(4,:),plotpar,sfac);
 eldisp2(ex5,ey5,Ed(5,:),plotpar,sfac);
 eldisp2(ex6,ey6,Ed(6,:),plotpar,sfac);
 eldisp2(ex7,ey7,Ed(7,:),plotpar,sfac);
 eldisp2(ex8,ey8,Ed(8,:),plotpar,sfac);
 eldisp2(ex9,ey9,Ed(9,:),plotpar,sfac);
 axis([-5 15 -1 15]); 
 pltscalb2(sfac,[1e-2 5 0]);
 axis([-5 15 -1 15]);
 title('displacements')

%----- Draw normal force diagram --------------------------------

 figure(2)
 plotpar=[2 1];
 sfac=scalfact2(ex3,ey3,es3(:,1),0.2);
 eldia2(ex1,ey1,es1(:,1),plotpar,sfac);
 eldia2(ex2,ey2,es2(:,1),plotpar,sfac);
 eldia2(ex3,ey3,es3(:,1),plotpar,sfac);
 eldia2(ex4,ey4,es4(:,1),plotpar,sfac);
 eldia2(ex5,ey5,es5(:,1),plotpar,sfac);
 eldia2(ex6,ey6,es6(:,1),plotpar,sfac);
 eldia2(ex7,ey7,es7(:,1),plotpar,sfac);
 eldia2(ex8,ey8,es8(:,1),plotpar,sfac);
 eldia2(ex9,ey9,es9(:,1),plotpar,sfac);
 axis([-5 15 -1 15]);
 pltscalb2(sfac,[3e4 5 0]);
 title('normal force')

%----- Draw shear force diagram ---------------------------------

 figure(3)
 plotpar=[2 1];
 sfac=scalfact2(ex5,ey5,es5(:,2),0.2);
 eldia2(ex1,ey1,es1(:,2),plotpar,sfac);
 eldia2(ex2,ey2,es2(:,2),plotpar,sfac);
 eldia2(ex3,ey3,es3(:,2),plotpar,sfac);
 eldia2(ex4,ey4,es4(:,2),plotpar,sfac);
 eldia2(ex5,ey5,es5(:,2),plotpar,sfac);
 eldia2(ex6,ey6,es6(:,2),plotpar,sfac);
 eldia2(ex7,ey7,es7(:,2),plotpar,sfac);
 eldia2(ex8,ey8,es8(:,2),plotpar,sfac);
 eldia2(ex9,ey9,es9(:,2),plotpar,sfac);
 axis([-5 15 -1 15]);
 pltscalb2(sfac,[3e4 5 0]);
 title('shear force') 

%----- Draw moment diagram --------------------------------------

 figure(4)
 plotpar=[2 1];
 sfac=scalfact2(ex4,ey4,es1(:,3),0.1);
 eldia2(ex1,ey1,es1(:,3),plotpar,sfac);
 eldia2(ex2,ey2,es2(:,3),plotpar,sfac);
 eldia2(ex3,ey3,es3(:,3),plotpar,sfac);
 eldia2(ex4,ey4,es4(:,3),plotpar,sfac);
 eldia2(ex5,ey5,es5(:,3),plotpar,sfac);
 eldia2(ex6,ey6,es6(:,3),plotpar,sfac);
 eldia2(ex7,ey7,es7(:,3),plotpar,sfac);
 eldia2(ex8,ey8,es8(:,3),plotpar,sfac);
 eldia2(ex9,ey9,es9(:,3),plotpar,sfac);
 axis([-5 15 -1 15]);
 pltscalb2(sfac,[10e4 5 0]);
 title('moment') 

%------------------------ end -----------------------------------
 echo off
jonaslindemann commented 2 years ago

You need to add an additional rotational degree of freedom at the ball joint to implement the ball joint.

        24 25
o---------o o---------o

You probably need to renumber the degrees of freedom.

jonaslindemann commented 2 years ago

No response. Closing.