% ------------------------------------------------------------- % Symbolic FEM ----------------------------------------------- % Djebar BAROUDI, PhD. % 27.3.2017 % ... to obtain FE, just write the correct integration scheme now to ... % ... integrate exactly. % --------------------------------------------------------------------- syms EK K syms EK_f EF_f EF syms SK SF syms SK_f SF_f syms SK_BC SK_BC_f SF_BC SF_BC_f syms K11 K12 K21 K22 syms EI L P P0 syms x syms h wk syms x_0 x_L x_R % xL := x_i-1, x0 = x_i, x_i+1 = xR syms L1 L2 L3 L4 % second order Lagrange polynoms L1,L2,L3:=L3.1, L3.2, L3.3 syms H1 H2 H3 H4 % cubic Hermite syms ksi % ksi = [-1 1], x = x1 + (1+ksi)2 * L, L = b-a syms a b syms x1 x2 % a = x1 , b = x2 % % external load syms q0 coeff % BCs syms v_BC syms v v_f v_f_disp v_f_rot txt = 'HUOM! Tässä symboli L := h' % % dx = L/2 * dksi --> df/dx = 2/L * df/dksi, d2f/dx2 = (2/L)**2 * d2f/dksi2 % % Hermite interpolation polynomes (or shape function) .. % ... within interval [a % b] % H1 --> f(a), H2 --> f'(a) , H3 --> f(b), H4 --> f'(b) H1(ksi, L) = ((1-ksi).^2).*(2+ksi)/4; H2(ksi, L) = L * (1-ksi.^2).*(1-ksi)/8; H3(ksi, L) = ((1+ksi).^2).*(2-ksi)/4; H4(ksi, L) = L *(-1+ksi.^2).*(1+ksi)/8; % 1st derivatives df/d_ksi of H_i, i = 1:3 d1_H1(ksi, L) = simplify( diff(H1, ksi) ); d1_H2(ksi, L) = simplify( diff(H2, ksi) ); d1_H3(ksi, L) = simplify( diff(H3, ksi) ); d1_H4(ksi, L) = simplify( diff(H4, ksi) ); % Huom! --> to obtain d_f/d_x = (2/L) * d_f/d_ksi % 2nd derivatives d2_f/d_ksi**2 of H_i, i = 1:3 d2_H1(ksi, L) = simplify( diff(d1_H1, ksi) ); d2_H2(ksi, L) = simplify( diff(d1_H2, ksi) ); d2_H3(ksi, L) = simplify( diff(d1_H3, ksi) ); d2_H4(ksi, L) = simplify( diff(d1_H4, ksi) ); % Huom! --> to obtain d2_f/d_x**2 = (2/L)^2 * d2_f/d_ksi**2 % % Element matrices ... % ------------------------ % Quadtrature weights wk = h * [1/2 1/2] ; ksi_k = [-1 +1]; % --------------------------------------------------- % Stifness matrix K_B (bending) % -------------------------------------------------- K11 = int( d2_H1(ksi, L) * d2_H1(ksi, L), ksi, [-1 +1] ); K12 = int( d2_H1(ksi, L) * d2_H2(ksi, L), ksi, [-1 +1] ); K13 = int( d2_H1(ksi, L) * d2_H3(ksi, L), ksi, [-1 +1] ); K14 = int( d2_H1(ksi, L) * d2_H4(ksi, L), ksi, [-1 +1] ); K21 = int(d2_H2(ksi, L) * d2_H1(ksi, L), ksi, [-1 +1]); K22 = int(d2_H2(ksi, L) * d2_H2(ksi, L), ksi, [-1 +1]); K23 = int(d2_H2(ksi, L) * d2_H3(ksi, L), ksi, [-1 +1]); K24 = int(d2_H2(ksi, L) * d2_H4(ksi, L), ksi, [-1 +1]); K31 = int(d2_H3(ksi, L) * d2_H1(ksi, L), ksi, [-1 +1]); K32 = int(d2_H3(ksi, L) * d2_H2(ksi, L), ksi, [-1 +1]); K33 = int(d2_H3(ksi, L) * d2_H3(ksi, L), ksi, [-1 +1]); K34 = int(d2_H3(ksi, L) * d2_H4(ksi, L), ksi, [-1 +1]); K41 = int(d2_H4(ksi, L) * d2_H1(ksi, L), ksi, [-1 +1]); K42 = int(d2_H4(ksi, L) * d2_H2(ksi, L), ksi, [-1 +1]); K43 = int(d2_H4(ksi, L) * d2_H3(ksi, L), ksi, [-1 +1]); K44 = int(d2_H4(ksi, L) * d2_H4(ksi, L), ksi, [-1 +1]); % ------------------------------------------------- % Element Stifness matrix K_B (bending) % int( f ).d_x = L/2 * int( f ). d_ksi sekä d2g/dx^2 = (2/L)^4 * d2g/d_ksi^2 coeff = ( ((2/L)^2) * ((2/L)^2) ) * L/2; EK_f(EI, L, h) = EI * coeff *[K11 K12 K13 K14 K21 K22 K23 K24 K31 K32 K33 K34 K41 K42 K43 K44]; EK = EI * coeff * [K11 K12 K13 K14 K21 K22 K23 K24 K31 K32 K33 K34 K41 K42 K43 K44]; % ---------------------------------------------------- % --------------------------------------------------- % Element loading vector d_ksi = 2/L dx % -------------------------------------------------- q(ksi, q0) = q0; % N/m f1 = L/2 * int( H1(ksi, L) * q(ksi, q0 ), ksi, [-1 +1]); f2 = L/2 * int( H2(ksi, L) * q(ksi, q0 ), ksi, [-1 +1]); f3 = L/2 * int( H3(ksi, L) * q(ksi, q0 ), ksi, [-1 +1]); f4 = L/2 * int( H4(ksi, L) * q(ksi, q0 ), ksi, [-1 +1]); % Element loading vector EF_f(h, q0) = [f1 f2 f3 f4]; EF = [f1 f2 f3 f4]; %% ---------------------------------------- %% Assembling global stifness matrix % --- IE = 1; %% NE = 4; NE = 8; NP = NE + 1; NDOFE = 4; % NDOFs / Element NCN = 2; % NDOFs / node %% NSDOF = 10; NSDOF = NCN*NP; %% NODOF = [1 2 3 4 %% 3 4 5 6 %% 5 6 7 8 %% 7 8 9 10]; NODOF = [1 2 3 4 3 4 5 6 5 6 7 8 7 8 9 10 9 10 11 12 11 12 13 14 13 14 15 16 15 16 17 18 ]; % --------------------------------------------- %% BEGIN: Assembly of the system matrices ---- % --------------------------------------------- SK = sym('SK', NSDOF); SK = sym(0*SK); SF = sym('SF', [1 NSDOF]); SF = sym(0*SF'); for IE = 1: NE for i = 1: NDOFE IS = NODOF(IE,i); SF(IS) = SF(IS) + EF(i); for j = 1 : NDOFE JS = NODOF(IE,j); SK(IS,JS) = SK(IS,JS) + EK(i, j); end end end % ----------------------------------- % Boundary conditions --- NBC = 2 %% NOBC = [1 9] ; % dofs with v(NOBC) = 0 NOBC = [1 17] ; % dofs with v(NOBC) = 0 v_BC(1) = 0; % the value of the DOF v_BC(2) = 0; % the value of the DOF v_BC = v_BC'; BC_GR = sym('BC_GR', NSDOF); BC_GR = sym(0*BC_GR) + 1; % --------------------- % BEGIN ---- BC % Exactly accounting for BCs in SK ----------- for IB = 1 : NBC IBC = NOBC(IB); BC_GR(IBC, 1: NSDOF) = 0; BC_GR(IBC, IBC) = 1; end SK_BC = BC_GR .* SK; % BCs in SK accounted ... now one should make SK symmetric SF_BC = SF; for IB = 1 : NBC IBC = NOBC(IB); SK_BC(IBC, IBC) = 1; % BCs account in the Stiffness matrix SF_BC(IBC) = v_BC(IB); % BCS account in the loading vector end SK_BC_f(EI, h, L) = SK_BC; % makes a function SF_BC_f(q0, h) = SF_BC; % ------------------------------------ % BCS account in the loading vector % ------------------------------------ % Note BENE : should make the global stiffness matrix symmetric (next time when I have time) % ------------------------- % Solution: % ---------------- v = inv(SK_BC) * SF_BC; txt = 'HUOM! Tässä symboli L := h , korvatiin L ---> L/NE' v = subs(v, L, L/NE); % because all above L: = h (the element length) % substitue to L --> L/NE, NE - element number v_f(EI, L, q0) = v; % solution v_f_disp(EI, L, q0) = v(1:2:NSDOF) v_f_rot(EI, L, q0) = v(2:2:NSDOF) v_tarkka_at_middle = 5/384 % * q0 *L^4 / EI % --------------------- % END ---- BC % Exactly accounting for BCs in SK ----------- % ------------------------------------------------- % Ploting shape functions - checking --------------------------------- % 26.3.2018 --- % --------------------------------------------- %% x0 = -1:0.01:1; x_00 = 0; x_LL = -1; x_RR = 1; %% L1_v = eval( L1(x0, x_00, x_LL, x_RR) )' figure fplot(H1(ksi, 1), [-1, 1],'r-') grid on hold on fplot(H2(ksi, 1), [-1, 1],'k-') fplot(H3(ksi, 1), [-1, 1],'b-') fplot(H4(ksi, 1), [-1, 1],'m-') xlabel('x') ylabel('Hermite interpolation polynomes H_i') legend('H_1','H_2','H_3','H_4') figure fplot(d1_H1(ksi, 1), [-1, 1],'r-') grid on hold on fplot(d1_H2(ksi, 1), [-1, 1],'k-') fplot(d1_H3(ksi, 1), [-1, 1],'b-') fplot(d1_H4(ksi, 1), [-1, 1],'m-') legend('dH_1','dH_2','dH_3','dH_4') figure fplot(d2_H1(ksi, 1), [-1, 1],'r-') grid on hold on fplot(d2_H2(ksi, 1), [-1, 1],'k-') fplot(d2_H3(ksi, 1), [-1, 1],'b-') fplot(d2_H4(ksi, 1), [-1, 1],'m-') legend('d2H_1','d2H_2','d2H_3','d2H_4') % -----------------