feat: Strain profile n my mz

[MestreCarlos]

A first draft for calculate the strain profile fon N+My+Mz forces:

• Exceptions are still to be launched when the stresses exceed the capacity of the section.
• The algorithm has room for optimisation because it works slowly
• Pending to create the test files

Example:

from shapely import Polygon

from structuralcodes import codes, materials
from structuralcodes.geometry import SurfaceGeometry
from structuralcodes.sections._generic import GenericSection
from structuralcodes.sections._reinforcement import add_reinforcement_line

# from structuralcodes.plots.section_plots import draw_section_response,draw_section,get_stress_point

### CREATE THE SECTION ########
codes.set_design_code(design_code='ec2_2004')
concrete = materials.concrete.create_concrete(fck=25)
reinforcemnet = materials.reinforcement.create_reinforcement(
    fyk=500,
    Es=200000,
    density=7850,
    ftk=500,
    epsuk=0.07,
)
poly = Polygon(((0, 0), (350, 0), (350, 500), (0, 500)))
geo = SurfaceGeometry(poly, concrete)
geo = add_reinforcement_line(
    geo, (50, 450), (300, 450), 12, reinforcemnet, n=3
)
geo = add_reinforcement_line(geo, (50, 50), (300, 50), 20, reinforcemnet, n=6)
sec = GenericSection(
    geo,
)
# sec.geometry = sec.geometry.translate(-sec.gross_properties.cy, -sec.gross_properties.cz)
sec.geometry = sec.geometry.translate(-175, -250)
### CREATE THE SECTION ########


### TEST calculate_strain_profile_biaxial  ########
n = -110 * 1e3
my = -31 * 1e6  # -150
mz = 65 * 1e6  # -100


res = sec.section_calculator.calculate_strain_profile(n, my, mz)
print(
    'eps',
    round(res[0] * 1e3, 2),
    '  chiy',
    round(res[1] * 1e6, 2),
    '  chiz',
    round(res[2] * 1e6, 2),
)
forces = sec.section_calculator.integrate_strain_profile(res)
print(
    'N',
    round(forces[0] / 1e3, 2),
    '  My',
    round(forces[1] / 1e6, 2),
    '  Mz',
    round(forces[2] / 1e6, 2),
)

[MestreCarlos]

I updated the branch to be reviewed. The method calculate_strain_profile works as follows:

1. First, check if the forces are beyond the section's capacity:

• Calculate the MyMz(N) diagram.
• Check if the point is outside the diagram using the nested function "check_points_in_2D_diagram." This function should be extracted for reuse. It is a generic function to evaluate whether a point is inside or outside a diagram.

2. We understand alpha to be the angle of the neutral axis and theta to be the angle of the resultant moments (My_ed, Mz_ed). I am looking for the angles (alpha_1 and alpha_2) that will define the range in which the algorithm will operate to determine the angle alpha that leads to My+Mz forces, integrated over the section, matching My_ed and Mz_ed.

3. I start the algorithm (bisection) by testing different values of alpha. For each alpha value, I calculate the moment-curvature diagram and look for the intersection point of the resultant moment to determine the total curvature. I decompose the total curvature into chi_y and chi_z using the alpha angle from the iteration. The moments My and Mz are obtained by interpolating chi_y and chi_z in the moment-curvature diagram.

Use example:

from shapely import Polygon
from structuralcodes import codes, materials
from structuralcodes.geometry import SurfaceGeometry
from structuralcodes.sections._generic import GenericSection
from structuralcodes.sections._reinforcement import add_reinforcement_line

# from structuralcodes.plots.section_plots import draw_section_response

# CREATE THE SECTION ########
codes.set_design_code(design_code='ec2_2004')
concrete = materials.concrete.create_concrete(fck=25)
reinforcemnet = materials.reinforcement.create_reinforcement(
    fyk=500,
    Es=200000,
    density=7850,
    ftk=500,
    epsuk=0.07,
)
poly = Polygon(((0, 0), (350, 0), (350, 500), (0, 500)))
geo = SurfaceGeometry(poly, concrete)
geo = add_reinforcement_line(
    geo, (50, 450), (300, 450), 12, reinforcemnet, n=3
)
geo = add_reinforcement_line(geo, (50, 50), (300, 50), 20, reinforcemnet, n=6)
sec = GenericSection(
    geo,
)
sec.geometry = sec.geometry.translate(-175, -250)
# CREATE THE SECTION ########


# TEST calculate_strain_profile_biaxial  ########
n = -260 * 1e3
my = 25 * 1e6  # -150
mz = 80 * 1e6  # -100

res = sec.section_calculator.calculate_strain_profile(n, my, mz)
print(
    'eps',
    round(res[0] * 1e3, 2),
    '  chiy',
    round(res[1] * 1e6, 2),
    '  chiz',
    round(res[2] * 1e6, 2),
)
forces = sec.section_calculator.integrate_strain_profile(res)
print(
    'N',
    round(forces[0] / 1e3, 2),
    '  My',
    round(forces[1] / 1e6, 2),
    '  Mz',
    round(forces[2] / 1e6, 2),
)

[talledodiego]

Thanks @MestreCarlos for the contribution! I will take a look to it in the coming days.

[talledodiego]

Thanks again @MestreCarlos for the contribution.

Seeking how to formulate a more robust and efficient algorithm I would propose the following idea based on Newton-Rhapson method (or similar versions like initial tangent or secant).

Determination of strain plain given external forces N, My, Mz

The vector containing the quantities $N$, $M_y$, $M_z$ is obtained by integrating the stresses over the cross-sectional area of a structural element.
They can be expressed as follows:

$$N = \int_A \sigma(y, z) \,dA$$ $$M_y = \int_A \sigma(y, z) \cdot z \,dA$$ $$M_z = \int_A \sigma(y, z) \cdot \left(-y\right) \,dA$$

or in vector form as:

$$\begin{bmatrix} N \\ M_y \\ M_z \end{bmatrix} = \int_A \begin{bmatrix} 1 \\\ z \\\ -y \end{bmatrix} \sigma\left(y, z\right) \, dA$$

The strain $\varepsilon(y,z)$ can be expressed as a function of the axial deformation and the curvatures as:

$$\varepsilon(y,z) = \varepsilon_0 + \kappa_y \cdot z - \kappa_z \cdot y$$

The stress $\sigma(y,z)$ can then be related to the strain $\varepsilon(y,z)$ through the constitutive law (in general non-linear), according to:

$$\sigma\left(y,z\right) = f\left[ \varepsilon(y,z) \right] = f \left( \varepsilon_0 + \kappa_y \cdot z - \kappa_z \cdot y \right)$$

Substituting this expression into the first equation the external forces can be written in matrix form as:

$$\begin{bmatrix} N \\ M_y \\ M_z \end{bmatrix} = \int_A \begin{bmatrix} 1 \\\ z \\\ -y \end{bmatrix} f \left( \begin{bmatrix} 1 & z & -y \end{bmatrix} \cdot \begin{bmatrix} \varepsilon_0 \\\ \kappa_y \\\ \kappa_z \end{bmatrix} \right) \, dA = \int_A \begin{bmatrix} 1 \\\ z \\\ -y \end{bmatrix} f \left( \begin{bmatrix} 1 & z & -y \end{bmatrix} \cdot \boldsymbol{\varepsilon} \right) \, dA$$

$\mathbf{x}$ where $\boldsymbol{\varepsilon}$ is the vector containing the strain plane coefficients. To solve for the unknowns coefficients given the external forces $\textbf{F}^{ext} = \left[N, \ M_y, \ M_z\right]^T$ we can use a numeric method like Newton-Raphson.

To this aim we can write the residual vector as:

$$\mathbf{R} \left( \mathbf{\varepsilon} \right) = \mathbf{F}^{ext} - \int_A \begin{bmatrix} 1 \\\ z \\\ -y \end{bmatrix} f \left( \begin{bmatrix} 1 & z & -y \end{bmatrix} \cdot \boldsymbol{\varepsilon} \right) \, dA$$

And linearizing it:

$$\mathbf{R} \left( \mathbf{\varepsilon}^{k+1} \right) = \mathbf{R} \left( \mathbf{\varepsilon}^k \right) + \mathbf{J} \left( \mathbf{\varepsilon}^k \right) \, \Delta\mathbf{\varepsilon}$$

It is therefore possible to solve for the current step:

$$\Delta\mathbf{\varepsilon} = -\mathbf{J} \left( \mathbf{\varepsilon}^k \right)^{-1} \, R \left( \mathbf{\varepsilon}^k \right)$$

writing so:

$$\mathbf{\varepsilon}^{k+1} = \mathbf{\varepsilon}^{k} + \Delta\mathbf{\varepsilon}$$

This iteration can be done until for instance the residual vector norm $\left | \left | R\left(\varepsilon^{k+1}\right)\right |\right |$ is below a selected tolerance.

The Jacobian matrix is basically the tangent stiffness of the section and is a 3x3 matrix given by:

$$\mathbf{J} = \int_A \frac{d\sigma}{d\mathbf{\varepsilon}} \begin{bmatrix} 1 \\\ z \\\ -y \end{bmatrix} \, dA$$

where:

$$\frac{d\sigma}{d\mathbf{\varepsilon}} = \frac{d\sigma}{d\varepsilon} \cdot \frac{d\varepsilon}{d\mathbf{\varepsilon}} = \frac{d\sigma}{d\varepsilon} \begin{bmatrix} 1 & z & -y \end{bmatrix}$$

Therefore the Jacobian matrix can be written as:

$$\mathbf{J} = \int_A \frac{d\sigma}{d\varepsilon} \begin{bmatrix} 1 \\\ z \\\ -y \end{bmatrix} \begin{bmatrix} 1 & z & -y \end{bmatrix} dA = \int_A \frac{d\sigma}{d\varepsilon} \begin{bmatrix} 1 & z & -y \\\ z & z^2 & -yz \\\ -y & -yz & y^2 \end{bmatrix} dA$$

The derivative $\frac{d\sigma}{d\varepsilon}$ represents the tangent modulus of the constitutive law.

Further remarks

Of course we could check first as you are proposing in your code if the point is inside or outside the domain before starting the iterative process. I would suggest for now leaving out this part letting the algorithm fail for the lack of convergence. In a following PR we can address this topic more consistently with the results objects.

Said that, are you willing to try exploring drafting this algorithm? Please let me know if I can assist anyhow in the process.

PS:
As a last comment, I have seen that you had to update the reinforced concrete property in CompundGeometry after rotation or translation: good that you found the bug (thanks!) since in my intended behavior it should have been unnecessary to manually update that. This was due to a small bug resolved not in the CompundGeometry but in the SurfaceGeometry in PR #203, merged already with dev.

[mortenengen]

Ref. our chat in another channel @DanielGMorenaFhecor, @talledodiego is happy to continue the implementation of this feature from where @MestreCarlos left it and along the lines sketched above. I am looking forward to seeing more progress on this one 😃

[talledodiego]

Ref. our chat in another channel @DanielGMorenaFhecor, @talledodiego is happy to continue the implementation of this feature from where @MestreCarlos left it and along the lines sketched above. I am looking forward to seeing more progress on this one 😃

As a Christmas present: done! ready for review!

[mortenengen]

Great job on this one @MestreCarlos and @talledodiego! Thanks for contributing!