1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Shear Flow Box Beam (Close to the answer)

  1. Apr 12, 2014 #1
    1. The problem statement, all variables and given/known data
    Untitled.png


    2. Relevant equations
    The moment and shear flow equations

    Solved with Matlab


    3. The attempt at a solution

    I am getting real close to the answer but something is off and I can't figure it out.

    Untitled1.png

    Untitled1.png


    The full code (Matlab)
    Code (Text):

    deg = pi/180;
    G   = 4.e6;

    %...Flange areas (in^2):
    Af = [1 0.5 0.5 1 1 0.5 0.5 1];
    nstringers = length(Af);

    %...Flange coordinates (in):
    z = [ 0 6 14 24 24 14 6 0];
    y = [ 0  0  0  0  10  10  10 10];

    %...Wall thicknesses (in):
    t = [0.04 0.04 0.04 0.05 0.04 0.04 0.04 0.05 0.025 0.025];
    npanels = length(t);

    %...Wall connectivity:
    node = [1 2
            2 3
            3 4
            4 5
            5 6
            6 7
            7 8
            8 1
            7 2
            6 3];

    %...Wall lengths:
    for i = 1:npanels
        s(i) = norm([y(node(i,2)) - y(node(i,1)) z(node(i,2)) - z(node(i,1))]);
    end

    %...Locate the centroid:
    zG = sum(z.*Af)/sum(Af);
    yG = sum(y.*Af)/sum(Af);

    %...Compute the centroidal moments of inertia:
    Iy  = sum(Af.*(z - zG).^2);
    Iz  = sum(Af.*(y - yG).^2);
    Iyz = sum(Af.*(y - yG).*(z - zG));

    %...Flange load gradients, Equation [4.8.2]:
    syms Vy Vz
    dPdx = 1/(Iy*Iz - Iyz^2)...
                           *( (Iy*Vy - Iyz*Vz)*(y - yG)...
                             +(Iz*Vz - Iyz*Vy)*(z - zG)).*Af;

    %...Stringer equilibrium at flanges 1 through 7:
    syms  q1  q2  q3  q4  q5  q6  q7  q8  q9  q10
    eq1 = q1 - q8 -       dPdx(1);
    eq2 = q2 - q1 -  q9 - dPdx(2);
    eq3 = q3 - q2 - q10 - dPdx(3);
    eq4 = q4 - q3 -       dPdx(4);
    eq5 = q5 - q4 -       dPdx(5);
    eq6 = q6 - q5 + q10 - dPdx(6);
    eq7 = q7 - q6 +  q9 - dPdx(7);

    %...Moment equivalence about flange 8, including dummy torque T:
    syms T
    A1 = 1/2*s(1)*s(8);
    A2 = 1/2*s(2)*s(8);
    A3 = 1/2*s(3)*s(8);
    A4 = 1/2*s(4)*(s(1) + s(2) + s(3));
    A9 = 1/2*s(9)*s(1);
    A10 = 1/2*s(10)*(s(1) + s(3));
    eq8 = 2*A1*q1 + 2*A2*q2 + 2*A3*q3 + 2*A4*q4 - 2*A9*q9 - 2*A10*q10 - ...
          (-Vy*14 - Vz*0) - T;

    %...Solve the above eight equations for the shear flows q1 through q8
    %   in terms of the redundants q9 & q10, the applied loads Vy & Vz and
    %   the dummy torque T:
    solution = solve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8, ...
                      q1, q2, q3, q4, q5, q6, q7, q8);

    %...Store the results (Equations [f] of the text) in the vector 'q':
    q = [solution.q1 solution.q2 solution.q3 solution.q4 ...
         solution.q5 solution.q6 solution.q7 solution.q8 ...
         q9 q10];

    %...Obtain the virtual shear flows from the actual ones by zeroing the
    %   applied loads and replacing the true redundant shear flows by
    %   their virtual counterparts and the dummy torque T by the virtual
    %   torque dT:
    syms dq9 dq10 dT
    dq = subs(q, {q9 q10 Vy Vz T}, {dq9 dq10 0 0 dT});

    %...Substitute the values of the external loads into the shear flow
    %   expressions:
     q = subs(q, {Vy Vz T}, {-8000 0 0});

    %...Calculate the internal complementary virtual work:
    syms L
    dWint = -L/G*sum(s./t.*q.*dq);

    %...The external complementary virtual work:
    syms theta
    dWext = theta*dT;

    %...PCVW:
    dW = dWint + dWext;

    %...This results in an equation of the form  c1*dq9 + c2*dq10 + c3*dT = 0.
    %   Isolate c1, c2 and c3:
    c1 = subs(dW, {dq9 dq10 dT}, {1 0 0});
    c2 = subs(dW, {dq9 dq10 dT}, {0 1 0});
    c3 = subs(dW, {dq9 dq10 dT}, {0 0 1});

    %...Solve these three equations for q9, q10 and theta:
    solution = solve(c1, c2, c3, q9, q10, theta);

    %...Substitute the shear flow results into the expressions for
    %   q1 through q8:
    q     = subs(q,{q9 q10},{solution.q9 solution.q10});

    theta = solution.theta;

    fprintf('\n----------------------------------------------------\n')
    disp(probTitle)
    fprintf('\n %s\n %s\n\n', programmer, date)
    disp(probStatement)

    fprintf('\n Flange   z(in)    y(in)   A(in^2)')
    for i = 1:8
        fprintf('\n %4.0f %8.2f %8.2f %8.2f',i,z(i),y(i),Af(i))
    end

    fprintf('\n\n Wall Thickness(in) Length(in)')
    for i = 1:10
        fprintf('\n%4.0f %8.2f       %8.4f',i,t(i),s(i))
    end

    fprintf('\n\n Centroid coordinates (in.)')
    fprintf('\n    zG = %g    yG = %g \n', zG, yG)

    fprintf('\n Moments of inertia (in.^4):')
    fprintf('\n    IGz = %g   IGy = %g    IGyz = %g \n', Iz, Iy, Iyz)

    fprintf('\n Shear loads:')
    fprintf('\n    Vy = -8000 lb     Vz = 0 lb \n')

    fprintf('\n Shear flows (lb/in):\n')
    for i = 1:10
        fprintf('    q%g = %10.4f\n', i, double(q(i)))
    end

    fprintf('\n Angle of twist:')
    fprintf('\n    theta/L = %g degrees/in', double(theta/deg/L))

    fprintf('\n----------------------------------------------------\n\n')
     

    Attached Files:

  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?



Similar Discussions: Shear Flow Box Beam (Close to the answer)
  1. Shear stress vector (Replies: 0)

Loading...