Shear Flow Box Beam (Close to the answer)

[zack7]

Homework Statement

Homework Equations

The moment and shear flow equations

Solved with Matlab

The Attempt at a Solution

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

The full code (Matlab)

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')

 

[KataKoniK]

I am not sure what the specific issue is with your code without seeing it, but here are a few general tips for troubleshooting:

1. Double check your equations: Make sure you have the correct equations for moment and shear flow, and that you are using the correct variables and units.

2. Check your input values: Make sure you are using the correct values for the flange areas, wall thicknesses, and wall lengths. Double check your calculations to ensure they are correct.

3. Use debugging tools: Many programming languages have debugging tools that allow you to step through your code and see the values of variables at each step. This can help you identify where the issue is occurring.

4. Print out intermediate values: You can use the "disp" function in Matlab to print out values at various points in your code. This can help you identify where the issue is occurring.

5. Use online resources: There are many online resources available for troubleshooting code, such as forums, tutorials, and step-by-step guides. You can also reach out to other scientists or programmers for help.