Disclaimer. The discussion reported below is nothing new and results are available in almost every book or textbook treating this topic. The main purpose, in this and in most of the other posts, is simply to clarify my ideas and to write down rather clear (to me, at least) and consistent notes, without any claim to correctness. On the other hand, in case you would find the content interesting and useful, I just kindly ask you to cite source and author. Of course, I am not a native English speaker, I just do my best to translate from Italian.

Introduction

When analyzing a reinforced concrete cross section subjected to simple bending (Fig. 1), under the usual reinforced concrete assumptions – e.g. plane sections remain plane – it often comes in handy, to calculate the bending moment M, simplifying the parabola-rectangular stress distribution through an equivalent plain rectangular (stress-block) distribution.

“Equivalent” means that the second distribution (the simpified one) originates the same resultant force and the same resultant moment compared to the first distribution. In other words, the integration of stresses acting on concrete in compression, and the location of the barycentre with respect to the neutral axis, must provide the same values in both cases.

In other words, the goal is the calculation of modulus |C| and location y_c of the resultant of compressive stresses, acting on concrete, via a couple of coefficients η and ξ that multiply either the compressed area (y_na·b) subjected to the design strength (f_cd), or the depth of the neutral axis (y_na), respectively:

\displaystyle \begin{aligned} &C = \eta \cdot \left( y_{na} \cdot b \cdot f_{cd} \right) &&\text{resultant compressive force} \\[3ex] &y_c = \xi \cdot y_{na} &&\text{depth of compressive force} \\[3ex] \end{aligned}

In Figure 1 symbols have the following meaning (C = T, being simple bending):

\begin{aligned} \small b &= \text{width of the section} \\[1ex] h &= \text{depth of the section} \\[1ex] d &= \text{effective depth of the section} \\[1ex] y_{na} &= \text{depth of the neutral axis} \\[1ex] \varepsilon_s &= \text{strain of steel in tension} \\[1ex] \varepsilon_c &= \text{strain of concrete at the compressed edge of the section} \\[1ex] \sigma_c &= \text{stress of concrete at the compressed edge of the section} \\[1ex] y_G &= \text{distance of compression centre of gravity from the neutral axis}\\[1ex] C &= \text{resultant compression force on concrete} \\[1ex] T &= \text{resultant tension force on steel} \\[1ex] \end{aligned} \\[3ex]

Stress-strain relationship for concrete

Concerning concrete, the assumed stress-strain (σ–ε) curve is parabola-rectangular (e.g. provided by Eurocode 2 for strength classes up to C50/60), shown in Figure 2. It comprises a parabola branch with downward concavity, passing through the origin and with vertex having abscissa ε₁ and ordinate f_cd, which corresponds to the design strength of concrete. For strain values ε > ε₁ and ε ≤ ε_u (ultimate strain), the stress value is constant and equal to f_cd (plateau). Incidentally, Eurocode 2 uses the symbols ε_c2 and ε_cu2 to denote the two key strain values.

The analytical expression of the parabola branch can be easily derived by forcing the passage through the origin (0;0) and through the point (ε₁; f_cd), and by imposing a null value for the derivative calculated in ε₁ as well:

\displaystyle \begin{aligned} &\sigma_p (\varepsilon) = -A\varepsilon^2 + B\varepsilon + C &\text{[assuming A > 0]} \end{aligned} \\[3ex] \begin{aligned} &I) &&\sigma_p(0) = 0 \quad \Rightarrow \quad C = 0 \\[3ex] &II) &&\sigma_p(\varepsilon_1) = f_{cd} \quad \Rightarrow \quad -A\varepsilon_1^2 + B\varepsilon_1 = f_{cd} \\[3ex] &III) && \left. \frac{d\sigma(\varepsilon)}{d\varepsilon} \right|_{\varepsilon_1} = 0 \quad \Rightarrow \quad -2A\varepsilon_1 + B = 0 \end{aligned} \\[3ex] B = 2A\varepsilon_1 \quad \Rightarrow \quad \text{[based on II]} -A\varepsilon_1^2 + 2A\varepsilon_1^2 = f_{cd} \\[3ex] A = \frac{f_{cd}}{\varepsilon_1^2}; \quad B = 2\frac{f_{cd}}{\varepsilon_1} \\[3ex]

Putting everything together, the following formula can be obtained to express the strain-stress relationship:

\color{royalblue} \\[3ex] \boxed{ \sigma(\varepsilon) = \begin{cases} \displaystyle -\frac{f_{cd}}{\varepsilon_1^2}\varepsilon^2 + 2\frac{f_{cd}}{\varepsilon_1}\varepsilon & 0 \leq \varepsilon \leq \varepsilon_1 \\ f_{cd} & \varepsilon_1 \lt \varepsilon \leq \varepsilon_u \\ \end{cases} } \\[3ex]

A)   Resultant of compression stresses C

The compression force per unit width (C/b) derives from the stress integration along concrete:

\displaystyle \frac{C}{b} = \int_{0}^{y_{na}} \sigma(y) dy \\[3ex]

Stress has been expressed as a function of strain ε, but in this case it should be written as a function of the distance y from the neutral axis. Since a linear strain distribution over the section was assumed (hypothesys of plane sections that remain plane), based on a simple proportion the following result can be obtained:

\displaystyle \frac{\varepsilon}{y} = \frac{\varepsilon_c}{y_{na}} \quad \Rightarrow \quad \varepsilon = \left( \frac{\varepsilon_c}{y_{na}} \right) y \\[3ex]

Now, the problem has to be divided into two separated cases, according to the upper strain value ε_c, which can either be (1) ε_c ≤ ε₁ or (2) ε₁ < ε_c ≤ ε_u. For the sake of simplicity, a new parameter is defined:

\displaystyle \alpha = \frac{\varepsilon_c}{\varepsilon_1} \\[3ex]

(1)   Strain ε_c ≤ ε₁, that is α ≤ 1

The integral to solve becomes:

\displaystyle \begin{aligned} \frac{C}{b} &= \int_{0}^{y_{na}} \left[ -\frac{f_{cd}}{\varepsilon_1^2} \left( \frac{\varepsilon_c}{y_{na}} \right)^2 y^2 + 2\frac{f_{cd}}{\varepsilon_1} \left( \frac{\varepsilon_c}{y_{na}} \right) y \right] dy = \\[3ex] &= \left[ -\frac{1}{3}\frac{f_{cd}}{\varepsilon_1^2} \left( \frac{\varepsilon_c}{y_{na}} \right)^2 y^3 + \frac{2}{2} \frac{f_{cd}}{\varepsilon_1} \frac{\varepsilon_c}{y_{na}} y^2 \right]_{y_{na}} = -\frac{1}{3}\frac{\varepsilon_c^2}{\varepsilon_1^2} f_{cd} y_{na} + \frac{\varepsilon_c}{\varepsilon_1} f_{cd} y_{na} = \\[3ex] &= f_{cd} y_{na} \left[ \alpha \left(1 - \frac{1}{3} \alpha \right) \right] \\[3ex] \end{aligned} \\ C_1 = \eta_1 \; y_{na} b \; f_{cd} \\[3ex] \boxed{ \eta_1 = \alpha \left(1 - \frac{1}{3} \alpha \right)} \\[3ex]

(2)   Strain ε₁ < ε_c ≤ ε_u, that is 1 < α ≤ ε_u/ε₁

The integral can be easily solved by dividing the contribution of the parabola branch, from 0 to y₁, and the contribution of the plateau, from y₁ to y_na. Therefore:

\displaystyle \begin{aligned} \frac{C}{b} &= \int_{0}^{y_1} \sigma(y) dy + f_{cd} \left( y_{na} - y_1 \right) = \\[3ex] &= \frac{2}{3} f_{cd} y_1 + y_{na} f_{cd} - y_1 f_{cd} = \\[3ex] &= f_{cd} y_{na} - \frac{1}{3} f_{cd} y_1 \\[3ex] \end{aligned} \\

Recalling the definition of α, it results:

\displaystyle y_1 = \frac{y_{na}}{\alpha} \\[3ex] \frac{C}{b} = \left( 1- \frac{1}{3\alpha} \right) f_{cd} y_{na} = \left( \frac{3\alpha -1}{3\alpha} \right) y_{na} f_{cd} \\[3ex] C_2 = \eta_2 \; y_{na} b \; f_{cd} \\[3ex] \boxed{ \eta_2 = \frac{3\alpha -1}{3\alpha}} \\[3ex]

B)   Depth of the centre of gravity of compression stresses y_c

The location of compression force C is generally expressed as distance from the compressed edge. The distance from the neutral axis, by definition, can be obtained from the ratio between first moment of area (static moment) S of the distribution of compressive stresses and area A of the distribution itself:

\displaystyle y_G \cdot A_{tot} = S_{tot} \quad \Rightarrow \quad y_G = \frac{S_{tot}}{A_{tot}} \\[3ex] y_G = \frac{\int_{0}^{y_{na}} \sigma(y) \cdot y \; dy}{\int_{0}^{y_{na}} \sigma(y) \; dy} \\[3ex]

Again, two cases can be distinguished, the first for ε_c ≤ ε₁, corresponding to α ≤ 1, and the second for ε₁ < ε_c ≤ ε_u, that is 1 < α ≤ ε_u/ε₁.

(1)   Strain ε_c ≤ ε₁, that is α ≤ 1

Calculations go on as follows:

\begin{aligned} y_G &= \frac{\displaystyle \int_{0}^{y_{na}} \left[ -\frac{f_{cd}}{\varepsilon_1^2} \left( \frac{\varepsilon_c}{y_{na}} y \right)^2 + 2 \frac{f_{cd}}{\varepsilon_1} \left( \frac{\varepsilon_c}{y_{na}} y \right) \right] \cdot y \; dy}{\displaystyle \int_{0}^{y_{na}} \left[ -\frac{f_{cd}}{\varepsilon_1^2} \left( \frac{\varepsilon_c}{y_{na}} y \right)^2 + 2 \frac{f_{cd}}{\varepsilon_1} \left( \frac{\varepsilon_c}{y_{na}} y \right) \right] \; dy} = \\[6ex] &= \frac{\displaystyle \int_{0}^{y_{na}} \left[ -\frac{f_{cd}}{\varepsilon_1^2} \left( \frac{\varepsilon_c^2}{y_{na}^2} y^3 \right) + 2 \frac{f_{cd}}{\varepsilon_1} \left( \frac{\varepsilon_c}{y_{na}} y^2 \right) \right] \; dy}{\displaystyle \int_{0}^{y_{na}} \left[ -\frac{f_{cd}}{\varepsilon_1^2} \left( \frac{\varepsilon_c}{y_{na}} y \right)^2 + 2 \frac{f_{cd}}{\varepsilon_1} \left( \frac{\varepsilon_c}{y_{na}} y \right) \right] \; dy} = \\[6ex] &= \frac{\displaystyle \left[ -\frac{1}{4} \frac{f_{cd}}{\varepsilon_1^2} \left( \frac{\varepsilon_c^2}{y_{na}^2} y^4 \right) + \frac{2}{3} \frac{f_{cd}}{\varepsilon_1} \left( \frac{\varepsilon_c}{y_{na}} y^3 \right) \right]_0^{y_{na}}}{\displaystyle \left[ -\frac{1}{3} \frac{f_{cd}}{\varepsilon_1^2} \left( \frac{\varepsilon_c^2}{y_{na}^2} y^3 \right) + \frac{2}{2} \frac{f_{cd}}{\varepsilon_1} \left( \frac{\varepsilon_c}{y_{na}} y^2 \right) \right]_0^{y_{na}}} = \\[6ex] &= \frac{\displaystyle \left[ -\frac{1}{4} \frac{\cancel{f_{cd}}}{\varepsilon_1^2} \left( \frac{\varepsilon_c^2}{\cancel{y_{na}^2}} y_{na}^{\cancel{4} \; 2} \right) + \frac{2}{3} \frac{\cancel{f_{cd}}}{\varepsilon_1} \left( \frac{\varepsilon_c}{\cancel{y_{na}}} y_{na}^{\cancel{3} \; 2} \right) \right]}{\displaystyle \left[ -\frac{1}{3} \frac{\cancel{f_{cd}}}{\varepsilon_1^2} \left( \frac{\varepsilon_c^2}{\cancel{y_{na}^2}} y_{na}^{\cancel{3} \; 2} \right) + \cancel{\frac{2}{2}} \frac{\cancel{f_{cd}}}{\varepsilon_1} \left( \frac{\varepsilon_c}{\cancel{y_{na}}} y_{na}^{\cancel{2}} \right) \right]} = \\[3ex] &= \frac{\displaystyle \left( \frac{2}{3} - \frac{1}{4}\alpha \right) \cancel{\alpha}}{\displaystyle \left( 1 - \frac{1}{3}\alpha \right) \cancel{\alpha}} y_{na} = \frac{8 - 3\alpha}{4 \left(3-\alpha \right)} y_{na} \\[6ex] \end{aligned}

Recalling the definition of α, it can be obtained:

\displaystyle y_G = \frac{\left( \frac{2}{3} - \frac{1}{4}\alpha \right) \cancel{\alpha}}{\left( 1 - \frac{1}{3}\alpha \right) \cancel{\alpha}} y_{na} = \frac{8 - 3\alpha}{4 \left(3-\alpha \right)} y_{na} \\[3ex]

Since the depth value we are looking for is given by the difference between y_na and the quantity deduced above (y_G), which expresses the distance of the centre of gravity from the neutral axis, the following expression can be obtained:

\displaystyle \boxed{ \xi_1 = \frac{4-\alpha}{4 \left( 3-\alpha \right)}} \\[3ex]

(2)   Strain ε₁ < ε_c ≤ ε_u, that is 1 < α ≤ ε_u/ε₁

Recalling the additive property of first moments of area, it follows that:

\begin{aligned} y_G &= \frac{S_{par} + S_{rec}}{A_{tot}} = \frac{\displaystyle A_{par} \cdot \left( 1 - \frac{3}{8} \right) y_1 + A_{rec} \cdot \frac{1}{2} \left( y_{na} + y_1 \right) }{\displaystyle A_{par} + A_{rec}} \\[3ex] &= \frac{\displaystyle \frac{2}{3} y_1 \cancel{f_{cd}} \left( \frac{5}{8} \right) y_1 + \left( y_{na} - y_1 \right) \cancel{f_{cd}} \frac{1}{2} \left( y_{na} + y_1 \right) }{\displaystyle \frac{2}{3} y_1 \cancel{f_{cd}} + \left( y_{na} - y_1 \right) \cancel{f_{cd}}} \\[3ex] &= \frac{\displaystyle \frac{1}{2} y_{na}^2 - \frac{1}{12} y_1^2}{\displaystyle y_{na} - \frac{1}{3} y_1} \\[3ex] &\text{poiché} \quad \frac{\varepsilon_c}{y_{na}} = \frac{\varepsilon_1}{y_1} \quad \Rightarrow \quad y_1 = \frac{\varepsilon_1}{\varepsilon_c} y_{na} = \frac{y_{na}}{\alpha} \\[3ex] y_G &= \frac{\displaystyle \frac{1}{2} y_{na}^{\cancel{2}} - \frac{1}{12} y_{na}^{\cancel{2}} \frac{1}{\alpha^2}}{\displaystyle \cancel{y_{na}} \; 1 - \frac{1}{3} \frac{\cancel{y_{na}} \; 1}{\alpha}} = \frac{\displaystyle \frac{1}{\cancel{12} \; 4 \alpha^{\cancel{2}}} \left( 6 \alpha^2 -1 \right)}{\displaystyle \frac{1}{\cancel{3\alpha}} \left( 3\alpha -1 \right) } y_{na} \\[3ex] &= \frac{6\alpha^2 -1}{4\alpha \left( 3\alpha - 1 \right)} y_{na} \\[3ex] \end{aligned}

As before, the depth of the centre of gravity of compression stresses is given by the difference between the depth of the neutral axis y_na and the distance of the centre of gravity from itself, that is y_G. It therefore results:

\displaystyle \boxed{ \xi_2 = \frac{6\alpha^2 - 4\alpha +1}{4\alpha \left( 3\alpha - 1 \right)}} \\[3ex]

Summary of the obtained results

In brief, the following formulas can be obtained to calculate the filling coefficient η and the depth coefficient ξ as a function of α, which expresses the ratio between strains ε_c and ε₁:

First case:   0 < α ≤ 1   (i.e.  0 < ε_c ≤ ε₁)

\displaystyle \\[3ex] \color{royalblue} \boxed{ \begin{aligned} \eta_1 &= \alpha \left(1 - \frac{1}{3} \alpha \right); &\xi_1 &= \frac{4-\alpha}{4 \left( 3-\alpha \right)} \end{aligned} } \\[3ex]

Second case:   1 < α ≤ ε_u/ε₁   (i.e.  ε₁ < ε_c ≤ ε_u)

\displaystyle \\[3ex] \color{royalblue} \boxed{ \begin{aligned} \eta_2 &= \frac{3\alpha -1}{3\alpha}; &\xi_2 &= \frac{6\alpha^2 - 4\alpha +1}{4\alpha \left( 3\alpha - 1 \right)} \end{aligned} } \\[3ex]

Results provided by formulas above, in case of ε₁ = 2·10⁻³ and ε_u = 3.5·10⁻³ (strength classes up to C50/60), are tabulated in Table 1 and graphically shown in Figure 3. It can be noted that the position of the centre of gravity does not vary much, from 0.333 (that is 1/3) for \alpha \to 0, up to 0.416. Conversely, as obvious, the filling coefficient ranges between 0 and 0.816 in correspondence to the ultimate strain value. For \alpha = 1 results are consistent with the well known properties of parabolas, i.e. \eta = 0.667 (that is 2/3) and \xi = 0.375 (that is centre of gravity at 3/8 \text{ and }5/8 from the edges).

A quick validation of analytical and numerical results can be found in the file below, elaborated by means of the excellent freeware software SMath Studio Desktop, which I already used in the past for educational purposes.

Table 1 Values of η e ξ as a function of α, in case of ε₁ = 2·10⁻³ and ε_u = 3.5·10⁻³ (strength classes up to C50/60)

α = εc / ε1	εc	η	ξ
0.10	0.2·10⁻³	0.097	0.336
0.20	0.4·10⁻³	0.187	0.339
0.30	0.6·10⁻³	0.270	0.343
0.40	0.8·10⁻³	0.347	0.346
0.50	1.0·10⁻³	0.417	0.350
0.60	1.2·10⁻³	0.480	0.354
0.70	1.4·10⁻³	0.537	0.359
0.80	1.6·10⁻³	0.587	0.364
0.90	1.8·10⁻³	0.630	0.369
1.00	2.0·10⁻³	0.667	0.375
1.10	2.2·10⁻³	0.697	0.381
1.20	2.4·10⁻³	0.722	0.388
1.30	2.6·10⁻³	0.744	0.394
1.40	2.8·10⁻³	0.762	0.400
1.50	3.0·10⁻³	0.778	0.405
1.60	3.2·10⁻³	0.792	0.410
1.70	3.4·10⁻³	0.804	0.414
1.75	3.5·10⁻³	0.810	0.416

Acknowledgements: images were created and elaborated with Inkscape and IrfanView. The chart in Figure 3 was obtained with LibreOffice Calc. Elaborations in the attached file were performed with SMath Studio Desktop. Finally, formulas have been inserted by means of the WordPress plugin \KaTeX (link).