Solving the Shifted Neutral Axis Problem

Question

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.

Answer 1

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]] /. γ -> 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"}
]}

Answer 2

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.