# Homework Help: Shear Flow Box Beam (Close to the answer)

1. Apr 12, 2014

### zack7

1. The problem statement, all variables and given/known data

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.

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

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

File size:
42.8 KB
Views:
223