## Balance Equations: Examples and Problems

### Example 1

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

where , , and 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:

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 metres at the top varying linearly to metres at the bottom. The beam is used to carry a concentrated load of , has a density , and is subjected to a gravity field with . Assuming that (a) the only nonzero stress component is , (b) is constant on every cross section perpendicular to , and (c) the undeformed and deformed coordinates coincide, find the distribution of the stress component  along the length of the beam.

##### Solution

The radius and the area vary with according to the equations:

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:

Simplifying and rearranging:

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

Even with the very simplified assumption that  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:

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 and is under a constant horizontal body force component . Assuming that (a) is the only nonzero stress component, (b) is only a function of , and (c) the undeformed and deformed coordinates coincide, find the distribution of 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, 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:

Assuming that is a function of only then:

Integrating the above equation yields:

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

Therefore, the stress distribution is:

Notice that if the end of the beam 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

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

then, find the values of , , , and 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:

Rearranging:

These equations have to be satisfied at every point inside the object, i.e., they have to be satisfied for all possible values of , , and . Therefore , , , and .

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:

The body forces vector applied is given by:

If the stress field is in equilibrium for all ,, and , determine the values of ,,,,.

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

The body forces vector applied is given by:

Show that the values of the constants ,,, and 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:

If the stress field is in equilibrium with a zero body forces vector field and the values of is equal to 5 when , show that the values of the constants , , and 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:

Find the equilibrium body forces vector and the constants and given that for , units and 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 .

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

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

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 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 and is used to carry a concentrated of 100N. at the right end. The beam is subjected to a distributed load of . Assuming that (a) the only nonzero stress component is , (b) is constant on every cross section perpendicular to , and (c) the undeformed and deformed coordinates coincide, find the distribution of the stress component along the length of the beam.

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

If the density in the reference configuration is given by the function units of density, find the expressions for and . For , draw the contour plots of , , and .