# Solving the Shifted Neutral Axis Problem

Problem Description:

In solving the "Shifted Neutral Axis" method for the stress distribution of a vertical steel vessel supported on a concrete foundation, we can derive the following expression:

$$\frac{\left(\gamma+\alpha\right)-\tan (\alpha)}{\left(\gamma+\alpha\right) \sec (\alpha )-\sin (\alpha )}=-\beta$$

where,

$$\alpha$$ is the angle of the "shifted" axis of rotation ranging from $$0$$ to $$\pi$$

$$\beta$$ is a ratio of dead to lateral load from $$0$$ to $$1$$. Real life values usually range from $$0.10$$ to $$0.50$$

$$\gamma$$ is a positive geometric and material constant. Real life values can range from $$0.05$$ to $$0.20$$

Essentially, this expression suggests that as the lateral load grows, the neutral axis shifts from a maximum value of $$\alpha=\pi$$ towards an asymptotic, smaller non-zero value.

Objective:

Usually, Engineers are given geometric and load constraints ($$\beta$$ and $$\gamma$$) and are told to find the system's response, $$\alpha$$.

The above expression, being implicit, leads itself to a simple solution via FindRoot[], but it would be ideal to solve for $$\alpha$$ as an explicit function of $$\beta$$ and $$\gamma$$, or at the very least come up with an approximated expression to calculate $$\alpha$$ directly.

I am not sure where to start solving this problem with Mathematica, but pointers would be appreciated.

Edit: The ultimate objective is to obtain an expression that reasonably approximates the shifted neutral axis angle--one that engineers can easily program into an Excel spreadsheet.

Generating a differential equation, you can calculate an interpolating function alphasol depending an beta and gamma.

Since i use version 8.0, i had to make an unusal construction for NDSolve. With higher versions use ParametricNDSolve with beta as parameter.

eq1[α_, β_, γ_] =
eq = ((γ + α) -
Tan[α])/((γ + α) Sec[α] -
Sin[α]) == -β

beta[α_] = -eq[] /. γ -> 1/10

invbeta[β_] = InverseFunction[beta][β]

deqγ =
D[eq1[α[β, γ], β, γ], γ] //
Simplify

αsol = α /.
First@NDSolve[{deqγ, α[β, 1/10] ==
invbeta[β]}, α, {β, .1, .5}, {γ, .05, \
.2}]


Get alpha-values directly, e.g. αsol[.1, .1] and compare with ContourPlot.

αsol[.1, .1]

{Manipulate[
ContourPlot[((γ + α) -
Tan[α])/((γ + α) Sec[α] -
Sin[α]) == -β, {γ, .05, 0.2}, {α, 0,
1.5}, AspectRatio -> 1, ImageSize -> 300], {β, .1, .5,
Appearance -> "Labeled"}],
Manipulate[
Plot[αsol[β, γ], {γ, .05, .2},
PlotRange -> {0, 1.5}, Frame -> True, AspectRatio -> 1,
ImageSize -> 300], {β, .1, .5, Appearance -> "Labeled"}
]}

• This is a neat way of getting plots from the interpolating function. The dream would be to generate an expression (analytically derived or fit) that other engineers can program into a spreadsheet. Jul 6, 2020 at 1:32

This is a transcendental equation and in general they don't have analytic solutions. There are special cases, of course, where a "trick" works. But in general no.

See https://en.wikipedia.org/wiki/Transcendental_equation for example.