Open Educational Resources

Balance Equations: Examples and Problems

Example 1

A stress field over a body that is in static equilibrium has the following form:

    \[\sigma =  \begin{pmatrix} 5x_1^2 + 3x_2 + x_3 & -x_2x_1 & 0 \\ -x_2x_1 & 5x_2 & 0 \\ 0 & 0 & 9x_3^2 \\ \end{pmatrix}\]

where x_1, x_2, and x_3 are the coordinates inside the body. Find the body forces vector field that is in equilibrium with this stress field.

Soultion

Since the given stress matrix is symmetric, it automatically satisfies the angular momentum balance equations. Using Eq. 2 above, the body forces vector field has the following form:

    \[\rho b =  \begin{pmatrix} -9x_1 \\ x_2 - 5 \\ -18x_3 \\ \end{pmatrix}\]

View Mathematica Code

x={x1,x2,x3};
s={{5*x1^2+3×2+x3,-x2*x1,0},{-x2*x1,5×2,0},{0,0,9×3^2}};
rhob1=-Sum[D[s[[i,1]],x[[i]]],{i,1,3}]
rhob2=-Sum[D[s[[i,2]],x[[i]]],{i,1,3}]
rhob3=-Sum[D[s[[i,3]],x[[i]]],{i,1,3}]

View Python Code

from sympy import Matrix,diff,symbols
import sympy as sp
sp.init_printing(use_latex="mathjax")
x1,x2,x3 = symbols("x1 x2 x3")
s = Matrix([[5*x1**2+3*x2+x3, -x2*x1, 0],[-x2*x1,5*x2,0],[0,0,9*x3**2]])
display("stress field over a body")
display("\u03C3 =", s)
x = Matrix([x1,x2,x3])
pb = Matrix([-sum([diff(s[i,j], x[i]) for i in range(3)]) for j in range(3)])
display("vector field")
display("pb =",pb) 

Example 2

The shown vertical beam has a varying circular cross sectional area with a radius r=2 metres at the top varying linearly to r=1 metres at the bottom. The beam is used to carry a concentrated load of P=50N, has a density \rho=20 kg/m^3, and is subjected to a gravity field with g = 10 m/sec^2. Assuming that (a) the only nonzero stress component is \sigma_{22}, (b) \sigma_{22} is constant on every cross section perpendicular to e_2, and (c) the undeformed and deformed coordinates coincide, find the distribution of the stress component \sigma_{22} along the length of the beam.

Solution

The radius r and the area A vary with X_2 according to the equations:

    \[r = (2-\frac{X_2}{L})\]

    \[A = \pi (2-\frac{X_2}{L})^2\]

A horizontal slice of the beam can be analyzed to find the equilibrium equation as follows:

Equating the sum of the vertical forces acting on the slice to zero yields:

    \[-\sigma_{22}A + \rho g A (dX_2) + \Big(\sigma_{22} + \frac{\partial \sigma_{22}}{\partial X_2}(dX_2)\Big)\Big(A + \frac{\partial A}{\partial X_2}(dX_2)\Big) = 0\]

Simplifying and rearranging:

    \[\frac{\partial \sigma_{22}}{\partial X_2}A + \sigma_{22}\frac{\partial A}{\partial X_2} + \rho g A = 0\]

A boundary condition for the stress is given at the free end as follows:

    \[@X_2 = L : \sigma_{22} = \frac{P}{A|_{X_2=L}} = \frac{P}{\pi}\]

Even with the very simplified assumption that \sigma_{22} is the only nonzero component of the stress matrix, the differential equation is relatively complicated. The command “DSolve” in Mathematica can be used to solve differential equations and was used to find the following distribution of the stress:

    \[\sigma_{22} = \frac{1}{3} (2\rho gL - \rho g X_2 + (\frac{L^2(3P-\rho g L \pi)}{\pi (X_2 - 2L)^2})) = \frac{1}{3} (2000 - \frac{23806.3}{(X_2 - 10)^2} - 200X_2) \hspace{1mm} N/m^2\]

View Mathematica Code

Clear[s,x,L,g]
A=Pi*(2-x/L)^2;
a=FullSimplify[DSolve[{s'[x]*A+D[A,x]*s[x]+ro*g*A==0,(s[L]==P/A/.x->L)},s[x],x]]
ss=s[x]/.a[[1]]
N[ss/.{L->5,ro->20,g->10,P->50}]

View Python Code

from sympy import diff,dsolve,simplify
import sympy as sp
x, rho, g, L, P = sp.symbols("x rho g L P")
f = sp.Function('f')
A = sp.pi*(2-x/L)**2
# sub in L into A
Al = A.subs({x:L})
display("area varying with X_2")
display("A =",A)
display("A x@L=",Al)
# partial derivative with respect to x, select f(L)=P/A 
sol = dsolve(f(x).diff(x)*A+A.diff(x)*f(x)+rho*g*A,f(x),ics={f(L):P/Al})
fsol=sol.rhs
display('solution f(x)=',fsol)
display("comparison with the Mathematica solution")
onethird = sp.Rational('1/3')
s = onethird*(2*rho*g*L-rho*g*x+L**2*(3*P-rho*g*L*sp.pi)/sp.pi/(x-2*L)**2)
display("Mathematica solution:",s)
dif = simplify(s-fsol)
display("difference:",dif)

Example 3

The shown horizontal beam has a density \rho=1 kg/m^3 and is under a constant horizontal body force component b_1=10 N/kg. Assuming that (a) \sigma_{11} is the only nonzero stress component, (b) \sigma_{11} is only a function of x_1, and (c) the undeformed and deformed coordinates coincide, find the distribution of \sigma_{11} along the beam.

Solution

It should be noted that the first two assumptions are very strong. For example, closer to the fixed boundary, there could be other nonzero stress components which would cause stress concentrations. Additionally, \sigma_{11} could be higher in the middle and lower closer to the boundaries. However, to easily achieve a closed form solution to the equilibrium equations, such assumptions are required. Otherwise, a numerical solution (for example, using finite element analysis) would be required to find all the components of the stress.

The equation of equilibrium in the horizontal direction is given by:

    \[\frac{\partial \sigma_{11}}{\partial x_1} + \rho b_1 = 0\]

Assuming that \sigma_{11} is a function of only x_1 then:

    \[\frac{d \sigma_{11}}{d x_1} + \rho b_1 = 0\]

Integrating the above equation yields:

    \[\sigma_{11} = -10x_1 + C\]

where C is a constant that can be obtained from the boundary condition given for the stress:

    \[@x_1 = 5: \sigma_{11} = 20 \Rightarrow -10(5) + C = 20 \Rightarrow C = 70N/m^2\]

Therefore, the stress distribution is:

    \[\sigma_{11} = -10x_1 + 70N/m^2\]

Notice that if the end of the beam x_1=5 meters was also fixed, then the given boundary conditions would not be sufficient to solve the differential equation of equilibrium and a constitutive law would be required to replace the stresses with expressions of the displacement in the differential equation of equilibrium.

View Python Code

from sympy import Eq,diff,integrate, dsolve
import sympy as sp
sp.init_printing(use_latex="mathjax")
rho, b1, x1, c, s = sp.symbols("rho b_1 x_1 C sigma")
f = sp.Function("f")
equation = f(x1).diff(x1)+rho*b1
display("equation of equilibrium =",Eq(equation,0))
#intEqu = integrate(equation, x1)
#intEqu = intEqu.subs({b1:10,rho:1})
sol = dsolve(equation.subs({b1:10,rho:1}),f(x1),ics={f(5):20})
display("solution for the stress f(x1) (N/m^2)",sol)

Example 4

The following is the stress field inside a body in static equilibrium

    \[\sigma =  \begin{pmatrix} \alpha x_1^2 + \delta x_1 + 3x_2 + x_3 & -\beta x_2 x_1 & 0 \\ -\beta x_2x_1 & 5x_2^2 & 0 \\ 0 & 0 & \gamma x_3^2 \\ \end{pmatrix}\]

If the body is subjected to the following body forces vector:

    \[\rho b =  \begin{pmatrix} 10x_1 \\ 10x_2 \\ 10x_3 \\ \end{pmatrix}\]

then, find the values of \alpha, \beta, \delta, and \gamma such that static equilibrium is achieved.

Solution

The first step is to check that the stress matrix is symmetric to satisfy balance of angular momentum. Indeed, the given stress matrix is symmetric.
The equations of static equilibrium are given as:

    \[\frac{\partial \sigma_{11}}{\partial x_1} + \frac{\partial \sigma_{21}}{\partial x_2} + \frac{\partial \sigma_{31}}{\partial x_3} + \rho b_1 = 2\alpha x_1 + \delta - \beta x_1 + 10x_1 = 0\]

    \[\frac{\partial \sigma_{12}}{\partial x_1} + \frac{\partial \sigma_{22}}{\partial x_2} + \frac{\partial \sigma_{32}}{\partial x_3} + \rho b_2 = -\beta x_2 + 10x_2 + 10x_2 = 0\]

    \[\frac{\partial \sigma_{13}}{\partial x_1} + \frac{\partial \sigma_{23}}{\partial x_2} + \frac{\partial \sigma_{33}}{\partial x_3} + \rho b_3 = 2\gamma x_3 + 10x_3 = 0\]

Rearranging:

    \[(2\alpha - \beta + 10)x_1 + \delta = 0\]

    \[(-\beta + 20)x_2 = 0\]

    \[(2\gamma + 10)x_3 = 0\]

These equations have to be satisfied at every point inside the object, i.e., they have to be satisfied for all possible values of x_1, x_2, and x_3. Therefore \alpha=5, \beta=20, \delta =0, and \gamma=-5.

View Mathematica Code

x= {x1, x2, x3};
s = {{alpha*x1^2 + delta*x1 + 3*x2 + x3, -beta*x2*x1, 0}, {-beta*x2*x1, 5 x2^2, 0}, {0, 0, gamma*x3^2}};
pb = {10 x1, 10 x2, 10 x3};
Print["Equilibrium Equations:"]
Eq1 = Sum[D[s[[i, 1]], x[[i]]], {i, 1, 3}] + pb[[1]] == 0
Eq2 = Sum[D[s[[i, 2]], x[[i]]], {i, 1, 3}] + pb[[2]] == 0
Eq3 = Sum[D[s[[i, 3]], x[[i]]], {i, 1, 3}] + pb[[3]] == 0
Print["Solving Equations:"]
ss1 = ForAll[{x1, x2, x3}, Eq1 && Eq2 && Eq3]
Solve[Resolve[ss1], {alpha, beta, delta, gamma}]

View Python Code

from sympy import diff,Eq,solve,Matrix
import sympy as sp
sp.init_printing(use_latex="mathjax")
x1,x2,x3 = sp.symbols("x1 x2 x3")
a,b,c,d = sp.symbols("\u03B1 \u03B2 \u03B3 \u03B4")
x = Matrix([x1,x2,x3])
s = Matrix([[a*x1**2+d*x1+3*x2+x3,-b*x2*x1,0],
            [-b*x2*x1, 5*x2**2,0],
            [0,0,c*x3**2]])
display("stress field inside a body in static equilibrium")
display("\u03C3 =", s)
pb = Matrix([x1, x2, x3])*10
display("body force vector")
display("pb =", pb)
# systems of equation:
# d/dx_1[x_11]+d/dx_2[x_21]+d/dx_3[x_31]+pb_1=0
# d/dx_1[x_12]+d/dx_2[x_22]+d/dx_3[x_32]+pb_2=0
# d/dx_1[x_13]+d/dx_2[x_23]+d/dx_3[x_33]+pb_3=0
equation = [sum([diff(s[i,j], x[i])for i in range(3)])+pb[j] for j in range(3)]
display("system of equations =",Eq(equation[0],0),Eq(equation[1],0),Eq(equation[2],0))
Eq0 = Eq(d,0)
Eq1 = Eq(equation[0].coeff(x1),0)
Eq2 = Eq(equation[1].coeff(x2),0)
Eq3 = Eq(equation[2].coeff(x3),0)
display("because these equations are valid for all values of x1, x2, and x3. They can be refined to =",Eq0,Eq1,Eq2,Eq3)
equations=[d,equation[0].coeff(x1),equation[1].coeff(x2),equation[2].coeff(x3)]
sol = solve(equations,[a,b,c,d])
display("Solution:")
display("\u03B1 =", sol[a])
display("\u03B2 =", sol[b])
display("\u03B3 =", sol[c])
display("\u03B4 =", sol[d])

Problems

  1. A stress field inside a continuum has the following form:

    \sigma = \begin{pmatrix} \alpha x_1^2 + x_2 & \beta(x_1-5)x_2 & 0 \\ \omega(x_1-5)x_2 & 2x_2^2 & 0 \\ 0 & 0 & \gamma x_3 \\ \end{pmatrix}

    The body forces vector applied is given by:

    \rho b = \begin{pmatrix} \delta \\ 0 \\ 0 \\ \end{pmatrix}

    If the stress field is in equilibrium for all x_1,x_2, and x_3, determine the values of \alpha,\beta,\gamma,\delta,\omega.

  2. A stress field inside a continuum has the following form:

    \sigma = \begin{pmatrix} \alpha x_1^2 + \delta x_1 + 3x_2 + x_3 & -\beta x_2x_1 & 0 \\ -\beta x_2x_1 & 5x_2^2 & 0 \\ 0 & 0 & \gamma x_3 \\ \end{pmatrix}

    The body forces vector applied is given by:

    \rho b = \begin{pmatrix} 10 \\ 0 \\ 0 \\ \end{pmatrix}

    Show that the values of the constants \alpha,\beta,\delta, and \gamma so that the stress field is in equilibrium with the body forces vector field are 5, 10, -10, and 0, respectively.

  3. A stress field inside a continuum has the following form:

    \sigma = \begin{pmatrix} \alpha x_1 (x_1-5) + \gamma & \beta x_1x_2(x_1-5)(x_2-5) & 0 \\ \beta x_1x_2(x_1-5)(x_2-5) & 0 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}

    If the stress field is in equilibrium with a zero body forces vector field and the values of \sigma_{11} is equal to 5 when x_1=0, show that the values of the constants \alpha, \beta, and \gamma are 0, 0, and 5, respectively.

  4. A two dimensional stress field on a rectangular plate of dimensions 1 and 5 units and 0.01 units of thickness has the following form:

    \sigma = \begin{pmatrix} 5x_1x_2 + B & A & 0 \\ A & 0 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}

    Find the equilibrium body forces vector and the constants A and B given that for x_1=0, \sigma_{11}=2 units and \sigma_{12}=2 units. Then, drawn the contour plot of the von Mises stress on the plate and identify the critical location (maximum von Mises stress) if the bottom left corner of the rectangular plate has coordinates (0,0).

  5. A stress field inside a continuum has the following form:

    \sigma = \begin{pmatrix} 2x_1(x_2-3) & 2x_2(x_1+3) & x_3^2 \\ 2x_2(x_1+3) & 2x_1x_2 & 4x_3 \\ x_3^2 & 4x_3 & 0 \\ \end{pmatrix}

    Assume the mass density of the continuum to be 2 units and that the body forces vector per unit mass is given by:

    \sigma = \begin{pmatrix} -x_1-x_2-x_3 \\ -x_1-x_2 \\ 0 \\ \end{pmatrix}

    Is the continuum in equilibrium? If not, what is the acceleration vector?

  6. The shown horizontal beam has a varying rectangular cross section with a constant thickness t of 1m. The height of the beam changes linearly from 1m. at the left end to 2m. at the right end. The beam has a density of 10kg/m^3 and is used to carry a concentrated of 100N. at the right end. The beam is subjected to a distributed load of q=100N/m. Assuming that (a) the only nonzero stress component is \sigma_{11}, (b)\sigma_{11} is constant on every cross section perpendicular to e_1, and (c) the undeformed and deformed coordinates coincide, find the distribution of the stress component \sigma_{11} along the length of the beam.

  7. The two dimensional displacement field of a 2 units by 3 units plate is given by:

        \[x_1 = X_1 + 0.01tX_2\]

        \[x_2 = X_2 + 0.02tX_1\]

    If the density in the reference configuration is given by the function \rho_0=1+X_1X_2 units of density, find the expressions for \rho_t(X,t) and \rho(x,t). For t=1, draw the contour plots of \det (F), \rho_t(X,1), and \rho(x,1).

Leave a Reply

Your email address will not be published.