Is it possible to add other Concrete Models in concreteproperties ?

[afugur]

Hello mr.@robbievanleeuwen,

I want to add Unconfined-Mander concrete model reference from Turkey Building Earthquake Code-2018. I tried to make this like EurocodeNonLinear class which it is in Concreteproperties library. Also I want this concrete model tension strength is zero. If I make this I take this error --> "brentq algorithm failed."

Is it possible to add this model in Concrete Properties ?

My Stress-Strain Profile Code:

class UnconfinedMander(ConcreteServiceProfile):

    strains: List[float] = field(init=False)
    stresses: List[float] = field(init=False)
    elastic_modulus: float
    ultimate_strain: float
    compressive_strength: float
    compressive_strain: float
    tensile_strength: float

    def __post_init__(
    self,
    ):

        self.strains = []
        self.stresses = []

        Esec = self.compressive_strength/0.002
        r = self.elastic_modulus/(self.elastic_modulus-Esec)

        fc_0035 = (self.compressive_strength * 2 * r) / (r-1+pow(2,r))

        zero_035 = np.linspace(0,0.0035,5,endpoint=False)
        r035_r05 = np.linspace(0.0035,0.005,5)

        self.strains.append(-0.005)
        self.stresses.append(0)

        for i in zero_035:

            self.strains.append(i)
            x = i/0.002
            fc = (self.compressive_strength * x * r) / (r-1+pow(x,r))
            self.stresses.append(fc)


        for i in r035_r05:
            self.strains.append(i)
            fc = ((0.005 - i) * fc_0035) / 0.0015
            self.stresses.append(fc)

Stress-strain profile graph is like this.

# Section
comp_str = 20 # Compressive Strength of Concrete
E = 5000 * math.sqrt(comp_str) # Elastic Modulus
# Geometry
b = 400 # Width
d = 400 # Height
cc = 30 # Concrete Cover
fi = 10 # Diameter of transverse
conc_unconfied = Concrete(
    name="Unconfined",
    density=2.4e-6,
    stress_strain_profile=ssp.UnconfinedMander(
        elastic_modulus=E,
        ultimate_strain=0.005,
        compressive_strength=comp_str,
        compressive_strain=0.002,
        tensile_strength=0,
    ),
    ultimate_stress_strain_profile=ssp.NonlinearUnconfined(
        elastic_modulus=E,
        ultimate_strain=0.005,
        compressive_strength=comp_str,
        compressive_strain=0.002,
        tensile_strength=0,
    ),
    alpha_squash=0.85,
    flexural_tensile_strength=3.5,
    colour="red",
)


steel_hd = SteelBar(
    name="Steel - Hardening",
    density=7.85e-6,
    stress_strain_profile=ssp.SteelHardening(
        yield_strength=500,
        elastic_modulus=200e3,
        fracture_strain=0.05,
        ultimate_strength=600,
    ),
    colour="grey",
)

col = rectangular_section(b, d,material=conc_unconfied)
col1 = rectangular_section(b-(cc+fi/2)*2, d-(cc+fi/2)*2,material=conc_unconfied).align_center(col)

col2 = (col - col1) + col1

geom = add_bar_rectangular_array(
    geometry=col2,
    area=201,
    material=steel_hd,
    n_x=2,
    x_s=304,
    n_y=2,
    y_s=304,
    anchor=(-152, -152),
    exterior_only=True,
)

geom.plot_geometry()

conc_sec = ConcreteSection(geom)
conc_sec.plot_section(title="")

res = conc_sec.moment_curvature_analysis()
res.plot_results()

scipy version ==> 1.6.2

[Agent6-6-6]

@afugur - FYI - Note I'm actually working at the moment on a specific implementation of this which relates to the New Zealand code implementation of the mander confined and unconfined concrete models.

I get the same warning when doing a moment curvature analysis (right at the start). Not looked into why yet, but I suspect it is possibly due to the brentq algorithm picking a value for the neutral axis that is beyond the failure strain of the concrete or steel possibly?

My curve for comparison, not 100% if the curve needs to return to zero or horizontal (was going to ask this question @robbievanleeuwen ?) as noted in the Eurocode curve a small horizontal return at both ends of the min/max strain regions had been added. (NZ implementation caps the max compressive strain at 0.004 which is why it does not return to spalling strain like your curve.

Note I tried your code @afugur and something must have changed in later versions of concrete-properties as the bars added were outside the section centred on 0,0, but section centroid was at 200,200 (noting you're still using alpha_squash which was removed). But after correcting that I see the same warning using your implementation.

[robbievanleeuwen]

@Agent6-6-6 I think this has been covered in previous comments e.g. in #47, but the horizontal stress is only for solver stability. You can never get a strain in the moment-curvature analysis above what get_ultimate_compressive_strain() gives you, but during the neutral axis finding algorithm you can have very small and large neutral axis depths, which result in large strains. So it is useful for the algorithm for the stresses to be in the ballpark of what it should be. If there was no horizontal stress you could conceivably get tensile stresses for large compressive strains after extrapolation, which can mess things up a bit!

[afugur]

@Agent6-6-6 firstly thanks for your answer. I use the v0.21 of concrete properties because of that I haven't removed the alpha squash. I agreed with this "Not looked into why yet, but I suspect it is possibly due to the brentq algorithm picking a value for the neutral axis that is beyond the failure strain of the concrete or steel possibly?" Maybe this depends on split_section_at_strains function.

def split_section_at_strains(
    concrete_geometries: List[Geometry],
    theta: float,
    point_na: Tuple[float],
    ultimate: bool,
    ultimate_strain: float = None,
    d_n: float = None,
    kappa: float = None,
) -> List[Geometry]:
    r"""Splits concrete geometries at discontinuities in its stress-strain profile.

    :param concrete_geometries: List of concrete geometries
    :type concrete_geometries: List[Geometry]
    :param float theta: Angle (in radians) the neutral axis makes with the horizontal
        axis (:math:`-\pi \leq \theta \leq \pi`)
    :param point_na: Point on the neutral axis
    :type point_na: Tuple[float]
    :param bool ultimate: If set to True, uses ultimate stress-strain profile
    :param float ultimate_strain: Concrete strain at failure
    :param float d_n: Depth of the neutral axis from the extreme compression fibre
    :param float kappa: Curvature

    :return: List of split geometries
    :rtype: List[:class:`sectionproperties.pre.geometry.Geometry`]
    """

    # create splits in concrete geometries at points in stress-strain profiles
    concrete_split_geoms = []

    for conc_geom in concrete_geometries:
        if ultimate:
            strains = (
                conc_geom.material.ultimate_stress_strain_profile.get_unique_strains()
            )
        else:
            strains = conc_geom.material.stress_strain_profile.get_unique_strains()

        # loop through intermediate points on stress-strain profile
        for idx, strain in enumerate(strains[1:-1]):
            # depth to point with `strain` from NA
            if ultimate:
                d = strain / ultimate_strain * d_n
            else:
                d = strain / kappa
            # convert depth to global coordinates
            dx, dy = global_coordinate(phi=theta * 180 / np.pi, x11=0, y22=d)

            # calculate location of point with `strain`
            pt = point_na[0] + dx, point_na[1] + dy

            # split concrete geometry (from bottom up)
            top_geoms, bot_geoms = split_section(
                geometry=conc_geom,
                point=pt,
                theta=theta,
            )

            # save bottom geoms
            concrete_split_geoms.extend(bot_geoms)

            # continue to split top geoms
            conc_geom = CompoundGeometry(geoms=top_geoms)

        # save final top geoms
        concrete_split_geoms.extend(top_geoms)

    return concrete_split_geoms

[Agent6-6-6]

@afugur, I figured out what the issue is:-

Once the brentq solver fails in concrete_section.moment_curvature_analysis, or more specifically when the service_normal_force_convergence method sets the _failure flag due to the defined limiting strains in either steel or concrete being exceeded. In this case the brentq solver returns a ValueError of f(a) and f(b) must have different signs error and it defaults to the except statement and raises the brentq algorithm failed. warning message as a result.

In the except statement it should first check for moment_curvature._failure = True, and only raise the ValueError if that is not true (as the error is then not the expected behaviour).

But effectively you can ignore the warning message, it's basically just saying the concrete or steel reached failure through failing to converge once the strains are exceeded.

I'll raise a pull request to address this.

[robbievanleeuwen]

@afugur thanks for bringing this up and thanks for solving @Agent6-6-6!

[robbievanleeuwen]

@Agent6-6-6, off-topic - wondering what part of NZ you are from? I'll be in Auckland next Feb for ~5 days and wondering if you have any must see recommendations?

[Agent6-6-6]

@robbievanleeuwen, yeah I'm based in Auckland. If you're sticking just in the greater Auckland region maybe look up some of the following and see if anything is to your liking:-

• Head up Skytower, or to One tree hill or Mt Eden to get a good overview of the city.
• If you're into walking there are plenty of walks from an hour to several hours all round Auckland, check out the Waitakere Ranges, or some of the coastal walks on east or west coast or some of the regional parks. https://www.aucklandnz.com/visit/discover/activities/walks-and-hikes?gclid=Cj0KCQjwhY-aBhCUARIsALNIC07fqxCr_ygF8vNCiLOJpRl-CLGoxlWkUEr-E4YkdsMUm3YZrUcQMnIaAlacEALw_wcB
• Head over to Waiheke Island via ferry for the day and do a winery tour.
• Take a ferry to Rangitoto and climb to the summit of the extinct volcano https://www.doc.govt.nz/parks-and-recreation/places-to-go/auckland/places/rangitoto-island/tracks/rangitoto-island-summit-walks/
• Check out some of the local beaches, there are some pretty rugged surf beaches on the west coast (Piha for example), but some of the inner city ones on north shore (Takapuna/Orewa) and eastern beachs (Mission bay) can be nice on a good day and have plenty of access to nice cafes for lunch/breakfast/dinner (should be nice weather in early February).
• Plenty of other touristy type places/activities within 1-3 hours travel if you have access to a car and want to get a bit further out of town to see more of the country. Rotorua (great for mountain biking if into that), Waitomo glo-worm Caves, Coromandel peninsula, etc.

Hope that helps!

[robbievanleeuwen]

@Agent6-6-6 much appreciated, thanks a ton 👍