diff --git a/index.ipynb b/index.ipynb
index 89aecdb41882b88a0f75052f16d271144f315a8a..386c97b050a9b443363e09e5a2c6976fa76ddc94 100644
--- a/index.ipynb
+++ b/index.ipynb
@@ -269,7 +269,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.9.1"
+ "version": "3.8.8"
},
"toc": {
"base_numbering": 1,
diff --git a/pull_out/2_1_1_PO_observation.ipynb b/pull_out/2_1_1_PO_observation.ipynb
index 909124aad96daaf887e88b83be5a8e2a0a31b250..2400abdedb18dad49e0641d105638e9c0ecc11c7 100644
--- a/pull_out/2_1_1_PO_observation.ipynb
+++ b/pull_out/2_1_1_PO_observation.ipynb
@@ -178,7 +178,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.9.1"
+ "version": "3.8.8"
},
"toc": {
"base_numbering": 1,
diff --git a/pull_out/2_2_1_PO_configuration_explorer.ipynb b/pull_out/2_2_1_PO_configuration_explorer.ipynb
index 7feacb69474c79dfc699008a656a832a6f910bdb..6380db65426b5bacc4f42402cd1c1c3f872c455c 100644
--- a/pull_out/2_2_1_PO_configuration_explorer.ipynb
+++ b/pull_out/2_2_1_PO_configuration_explorer.ipynb
@@ -280,7 +280,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.9.1"
+ "version": "3.8.8"
},
"toc": {
"base_numbering": 1,
diff --git a/pull_out/fragmentation.ipynb b/pull_out/fragmentation.ipynb
index d35f1d7d03f840cb86b246ba2b951cccf42f730a..713e7ffe99814b52c16b2ed6b2234ce6ee640837 100644
--- a/pull_out/fragmentation.ipynb
+++ b/pull_out/fragmentation.ipynb
@@ -47,14 +47,14 @@
},
{
"cell_type": "code",
- "execution_count": 7,
+ "execution_count": 4,
"id": "approved-pressing",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
- "model_id": "e338f59bf98748459c9275d04c2a6922",
+ "model_id": "42a4e4377cd5488f8e840614099ea927",
"version_major": 2,
"version_minor": 0
},
@@ -83,7 +83,7 @@
},
{
"cell_type": "code",
- "execution_count": 16,
+ "execution_count": 2,
"id": "olympic-validity",
"metadata": {},
"outputs": [
@@ -96,7 +96,7 @@
"E_\\mathrm{m}*Piecewise((P/(A_\\mathrm{f}*E_\\mathrm{f} + A_\\mathrm{m}*E_\\mathrm{m}), x <= -A_\\mathrm{m}*E_\\mathrm{m}*P/(A_\\mathrm{f}*E_\\mathrm{f}*\\bar{\\tau}*p + A_\\mathrm{m}*E_\\mathrm{m}*\\bar{\\tau}*p)), (-\\bar{\\tau}*p*x/(A_\\mathrm{m}*E_\\mathrm{m}), True))"
]
},
- "execution_count": 16,
+ "execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
@@ -238,7 +238,7 @@
},
{
"cell_type": "code",
- "execution_count": 9,
+ "execution_count": 3,
"id": "enabling-engine",
"metadata": {},
"outputs": [],
@@ -252,7 +252,7 @@
},
{
"cell_type": "markdown",
- "id": "interior-recycling",
+ "id": "0f66e3d8-e60d-4ff1-9a43-32ed2a486d7f",
"metadata": {},
"source": [
"The characteristic points of the ACK model deliver the values "
@@ -260,7 +260,7 @@
},
{
"cell_type": "code",
- "execution_count": 10,
+ "execution_count": 9,
"id": "sunset-trading",
"metadata": {},
"outputs": [],
@@ -292,7 +292,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
- "model_id": "6991363d3cc54733a4bf1d2e23c5cfe0",
+ "model_id": "ed30c8a8b1094c36b36199f212a8b434",
"version_major": 2,
"version_minor": 0
},
@@ -327,7 +327,7 @@
},
{
"cell_type": "code",
- "execution_count": 12,
+ "execution_count": 2,
"id": "boxed-scott",
"metadata": {
"tags": []
@@ -340,7 +340,7 @@
},
{
"cell_type": "code",
- "execution_count": 13,
+ "execution_count": 3,
"id": "regulated-republic",
"metadata": {
"tags": []
@@ -352,7 +352,7 @@
},
{
"cell_type": "code",
- "execution_count": 14,
+ "execution_count": 4,
"id": "subsequent-coach",
"metadata": {
"tags": []
@@ -361,7 +361,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
- "model_id": "eccd60851ba7434cadb9ad344207ccec",
+ "model_id": "63ba009ff6ac49c080b258484966fe0a",
"version_major": 2,
"version_minor": 0
},
@@ -393,10 +393,81 @@
"- Given the cross sectional areas $A_\\mathrm{m}, A_\\mathrm{f}$, concrete and reinforcement stiffness $E_\\mathrm{m}, E_\\mathrm{f}$ and strength $\\sigma_\\mathrm{mu}, \\sigma_\\mathrm{fu}$, reinforcement ratio $V_\\mathrm{f}$, bond intensity $T$, calculate the average crack width at failure of a tensile specimen."
]
},
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "26acec03-96b0-4f31-b890-4e6727d40060",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "VBox(children=(HBox(children=(VBox(children=(Tree(layout=Layout(align_items='stretch', border='solid 1px black…"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "ph.interact()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "694c155c-bc53-4750-82ef-a9b0b9f0a65d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%matplotlib widget\n",
+ "import pmcm.pmcm2"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "0d586905-6bf3-4f21-8a29-81973da643ff",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pm = pmcm.pmcm2.PMCM()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "11339d52-6dc2-49e0-9c60-7bf01b8d2b4b",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "137932dce7424c4091100ca25ef81d0b",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "VBox(children=(HBox(children=(VBox(children=(Tree(layout=Layout(align_items='stretch', border='solid 1px black…"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "pm.interact()"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
- "id": "second-entry",
+ "id": "be17e48c-81b5-493b-8e14-6f7631019f82",
"metadata": {},
"outputs": [],
"source": []
diff --git a/pull_out/pmcm/__init__.py b/pull_out/pmcm/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/pull_out/pmcm/api.py b/pull_out/pmcm/api.py
new file mode 100644
index 0000000000000000000000000000000000000000..c183405df21793738ccf3d0088b5f01d054b588a
--- /dev/null
+++ b/pull_out/pmcm/api.py
@@ -0,0 +1 @@
+from .pmcm import PMCM
\ No newline at end of file
diff --git a/pull_out/pmcm/pmcm.py b/pull_out/pmcm/pmcm.py
new file mode 100644
index 0000000000000000000000000000000000000000..be5586caf4e71d372f3efdf6cfd7b08e17a87cab
--- /dev/null
+++ b/pull_out/pmcm/pmcm.py
@@ -0,0 +1,276 @@
+"""
+Probabilistic multiple cracking model
+"""
+from typing import List, Any, Union
+import bmcs_utils.api as bu
+
+import traits.api as tr
+import numpy as np
+import warnings
+
+from bmcs_utils.trait_types import Float
+from scipy.optimize import newton
+
+warnings.filterwarnings("error", category=RuntimeWarning)
+
+
+class CrackBridgeModel(bu.Model):
+ """
+ Record of all material parameters of the composite. The model components
+ (PullOutModel, CrackBridgeRespSurf, PMCM) are all linked to the database record
+ and access the parameters they require. Some parameters are shared across all
+ three components (Em, Ef, vf), some are specific to a particular type of the
+ PulloutModel.
+ """
+
+class CBMConstantBond(CrackBridgeModel):
+ """
+ Return the matrix stress profile of a crack bridge for a given control slip
+ at the loaded end
+ """
+ Em = bu.Float(28000)
+ Ef = bu.Float(180000)
+ vf = bu.Float(0.01)
+ T = bu.Float(8)
+
+ ipw_view = bu.View(
+ bu.Item('Em'),
+ bu.Item('Ef'),
+ bu.Item('vf'),
+ bu.Item('T'),
+ )
+
+ @property
+ def Ec(self):
+ return self.Em * (1 - self.vf) + self.Ef * self.vf # [MPa] mixture rule
+
+class CrackBridgeRespSurface(bu.Model):
+ """
+ Crack bridge response surface that returns the values of matrix stress
+ along ahead of a crack and crack opening for a specified remote stress
+ and boundary conditions.
+ """
+ cb = bu.Instance(CrackBridgeModel)
+
+ def get_sig_m(self, z, sig_c):
+ """Get the profile of matrix stress along the specimen
+ :param z: np.ndarray
+ :type sig_c: float
+ """
+ cb = self.cb
+ sig_m = np.minimum(z * cb.T * cb.vf / (1 - cb.vf), cb.Em * sig_c /
+ (cb.vf * cb.Ef + (1 - cb.vf) * cb.Em))
+ return sig_m
+
+ def get_eps_f(self, z, sig_c):
+ cb = self.cb
+ sig_m = self.get_sig_m(sig_c, z )
+ eps_f = (sig_c - sig_m * (1 - cb.vf)) / cb.vf / cb.Ef
+ return eps_f
+
+ sig_c_slider: float = bu.Float(1.0, BC=True)
+
+ ipw_view = bu.View(
+ bu.Item('sig_c_slider')
+ )
+
+ @staticmethod
+ def subplots(fig):
+ ax = fig.subplots(1,1)
+ ax1 = ax.twinx()
+ return ax, ax1
+
+ def update_plot(self, axes):
+ ax, ax1 = axes
+ x_range = np.linspace(-1000,1000,1000)
+ z_range = np.abs(x_range)
+
+ sig_m_range = self.get_sig_m(self.sig_c_slider, z_range)
+ eps_f_range = self.get_eps_f(self.sig_c_slider, z_range)
+ sig_max = np.max(sig_m_range)
+ eps_max = np.max(eps_f_range)
+
+ ax.plot(x_range, sig_m_range, color='black')
+ ax.fill_between(x_range, sig_m_range, color='gray', alpha=0.1)
+ ax.set_ylim(ymin=-0.03*sig_max)
+ ax1.plot(x_range, eps_f_range, color='blue')
+ ax1.fill_between(x_range, eps_f_range, color='blue', alpha=0.1)
+ ax1.set_ylim(ymin=-0.03*eps_max)
+
+class PMCMHist(bu.Model):
+ pmcm = tr.WeakRef
+
+ K = bu.Int(0)
+ t = bu.Float(0)
+ K_max = tr.Property(depends_on='state_changed')
+ @tr.cached_property
+ def _get_K_max(self):
+ sig_c_K, eps_c_K, sig_mu_x, x, CS, sig_m_x_K, sig_m_x_K1 = self.pmcm.cracking_history
+ K_max = len(sig_c_K)
+ # ceil the index of current crack
+ self.K = np.min([self.K, K_max])
+ return K_max
+
+ ipw_view = bu.View(
+ bu.Item('K', latex=r'\mathrm{state}', readonly=True),
+ time_editor=bu.HistoryEditor(
+ var='t',
+ var_max='K_max'
+ )
+ )
+
+ @staticmethod
+ def subplots(fig):
+ ax1, ax2 = fig.subplots(1,2)
+ ax11 = ax1.twinx()
+ return ax1, ax11, ax2
+
+ def update_plot(self, axes):
+ ax, ax_cs, ax_sig_x = axes
+ cr = int(self.t * (self.K_max-2))
+ self.K = cr
+ self.pmcm.plot(axes)
+ sig_c_K, eps_c_K, sig_mu_x, x, CS, sig_m_x_K, sig_m_x_K1 = self.pmcm.cracking_history
+ ax_sig_x.plot(x,sig_m_x_K[cr])
+ ax_sig_x.plot(x,sig_m_x_K1[cr], linestyle='dashed')
+ ax.plot(eps_c_K[cr],sig_c_K[cr],color='magenta',marker='o')
+
+
+class PMCM(bu.Model):
+ name = "PMCM"
+ """
+ Implement the global crack tracing algorithm based on a crack bridge response surface
+ """
+
+ history = bu.Instance(PMCMHist)
+ def _history_default(self):
+ return PMCMHist(pmcm=self)
+
+ cb = bu.Instance(CrackBridgeModel)
+ def _cb_default(self):
+ return CBMConstantBond()
+
+ cb_rs = tr.Property#(depends_on="state_changed")
+ @tr.cached_property
+ def _get_cb_rs(self):
+ return CrackBridgeRespSurface(cb=self.cb)
+
+ tree = ['cb', 'cb_rs', 'history']
+
+ n_x = bu.Int(5000, ALG=True)
+ L_x = bu.Float(500, GEO=True)
+ sig_cu = bu.Float(20, MAT=True)
+ sig_mu = bu.Float(10, MAT=True)
+ m = bu.Float(4, MAT=True)
+
+ ipw_view = bu.View(
+ bu.Item('n_x'),
+ bu.Item('L_x'),
+ bu.Item('sig_cu'),
+ bu.Item('sig_mu'),
+ bu.Item('m'),
+ )
+
+ def get_z_x(self, x, XK): # distance to the closest crack (*\label{get_z_x}*)
+ """Specimen discretization
+ """
+ z_grid = np.abs(x[:, np.newaxis] - np.array(XK)[np.newaxis, :])
+ return np.amin(z_grid, axis=1)
+
+ def get_sig_c_z(self, sig_mu, z, sig_c_pre):
+ """
+ :param sig_c_pre:
+ :type sig_mu: float
+ """
+ # crack initiating load at a material element
+ #print('sig_mu', sig_mu)
+ fun = lambda sig_c: sig_mu - self.cb_rs.get_sig_m(sig_c, z)
+ try: # search for the local crack load level
+ sig_c = newton(fun, sig_c_pre)
+ #print('sig_c', sig_c)
+ return sig_c
+ except (RuntimeWarning, RuntimeError):
+ # solution not found (shielded zone) return the ultimate composite strength
+ return self.sig_cu
+
+ def get_sig_c_K(self, z_x, x, sig_c_pre, sig_mu_x):
+ # crack initiating loads over the whole specimen
+ get_sig_c_x = np.vectorize(self.get_sig_c_z)
+ sig_c_x = get_sig_c_x(sig_mu_x, z_x, sig_c_pre)
+ #print('sig_c_x', z_x, x, sig_c_pre, sig_mu_x)
+ #print('sig_c_x', sig_c_x)
+ y_idx = np.argmin(sig_c_x)
+ return sig_c_x[y_idx], x[y_idx]
+
+ cracking_history = tr.Property(depends_on='state_changed')
+ @tr.cached_property
+ def _get_cracking_history(self):
+ cb = self.cb
+ L_x, n_x, sig_mu, sig_cu, m = self.L_x, self.n_x, self.sig_mu, self.sig_cu, self.m
+ x = np.linspace(0, L_x, n_x) # specimen discretization
+ sig_mu_x: np.ndarray[np.float_] = sig_mu * np.random.weibull(
+ m, size=n_x) # matrix strength
+
+ XK: List[float] = [] # recording the crack postions
+ sig_c_K: List[float] = [0.] # recording the crack initiating loads
+ eps_c_K: List[float] = [0.] # recording the composite strains
+ CS: List[float] = [L_x, L_x / 2] # initial crack spacing
+ sig_m_x_K1: List[float] = [np.zeros_like(x)] # stress profiles at crack states
+ sig_m_x_K: List[float] = [np.zeros_like(x)] # stress profiles after crack states
+
+
+ Ec: float = cb.Ec
+ Em: float = cb.Em
+ idx_0 = np.argmin(sig_mu_x)
+ XK.append(x[idx_0]) # position of the first crack
+ sig_c_0 = sig_mu_x[idx_0] * Ec / Em
+ sig_m_x_K1.append(np.ones_like(x)*sig_mu_x[idx_0]) # matrix stress
+ #print('sig_c_0', sig_c_0)
+ sig_c_K.append(sig_c_0)
+ eps_c_K.append(sig_mu_x[idx_0] / Em)
+
+ while True:
+ z_x = self.get_z_x(x, XK) # distances to the nearest crack
+ sig_m_x_K.append(self.cb_rs.get_sig_m(sig_c_K[-1], z_x)) # matrix stress
+ sig_c_k, y_i = self.get_sig_c_K(z_x, x, sig_c_K[-1], sig_mu_x) # identify next crack
+ sig_m_x_K1.append(self.cb_rs.get_sig_m(sig_c_k, z_x)) # matrix stress
+ if sig_c_k == sig_cu:
+ break
+ XK.append(y_i) # record crack position
+ sig_c_K.append(sig_c_k) # corresponding composite stress
+ eps_c_K.append( # composite strain - integrate the strain field
+ np.trapz(self.cb_rs.get_eps_f(sig_c_k, self.get_z_x(x, XK)), x) / np.amax(x))
+ XK_arr = np.hstack([[0], np.sort(np.array(XK)), [L_x]])
+ CS.append(np.average(XK_arr[1:] - XK_arr[:-1])) # crack spacing
+
+ sig_c_K.append(sig_cu) # the ultimate state
+ eps_c_K.append(np.trapz(self.cb_rs.get_eps_f(sig_cu, self.get_z_x(x, XK) ), x) / np.amax(x))
+ CS.append(CS[-1])
+ return np.array(sig_c_K), np.array(eps_c_K), sig_mu_x, x, np.array(CS), np.array(sig_m_x_K), np.array(sig_m_x_K1)
+
+ @staticmethod
+ def subplots(fig):
+ ax1, ax2 = fig.subplots(1,2)
+ ax11 = ax1.twinx()
+ return ax1, ax11, ax2
+
+ def plot(self, axes):
+ ax, ax_cs, ax_sig_x = axes
+ sig_c_K, eps_c_K, sig_mu_x, x, CS, sig_m_x_K, sig_m_x_K1 = self.cracking_history
+ n_c = len(eps_c_K) - 2 # numer of cracks
+ ax.plot(eps_c_K, sig_c_K, marker='o', label='%d cracks:' % n_c)
+ ax.set_xlabel(r'$\varepsilon_\mathrm{c}$ [-]');
+ ax.set_ylabel(r'$\sigma_\mathrm{c}$ [MPa]')
+ ax_sig_x.plot(x, sig_mu_x, color='orange')
+ ax_sig_x.fill_between(x, sig_mu_x, 0, color='orange', alpha=0.1)
+ ax_sig_x.set_xlabel(r'$x$ [mm]');
+ ax_sig_x.set_ylabel(r'$\sigma$ [MPa]')
+ ax.legend()
+ eps_c_KK = np.array([eps_c_K[:-1], eps_c_K[1:]]).T.flatten()
+ CS_KK = np.array([CS[:-1], CS[:-1]]).T.flatten()
+ ax_cs.plot(eps_c_KK, CS_KK, color='gray')
+ ax_cs.fill_between(eps_c_KK, CS_KK, color='gray', alpha=0.2)
+ ax_cs.set_ylabel(r'$\ell_\mathrm{cs}$ [mm]');
+
+ def update_plot(self, axes):
+ self.plot(axes)
diff --git a/pull_out/pmcm/pmcm2.py b/pull_out/pmcm/pmcm2.py
new file mode 100644
index 0000000000000000000000000000000000000000..6c34898bde6d54d0f1397d91badec5dea1bc47d4
--- /dev/null
+++ b/pull_out/pmcm/pmcm2.py
@@ -0,0 +1,344 @@
+# %% md
+
+import numpy as np
+from scipy.optimize import newton
+import matplotlib.pylab as plt
+import traits.api as tr
+import bmcs_utils.api as bu
+
+# %%
+import warnings # (*\label{error1}*)
+
+warnings.filterwarnings("error", category=RuntimeWarning) # (*\label{error2}*)
+
+
+class CrackBridgeModel(bu.Model):
+ """
+ Record of all material parameters of the composite. The model components
+ (PullOutModel, CrackBridgeRespSurf, PMCM) are all linked to the database record
+ and access the parameters they require. Some parameters are shared across all
+ three components (Em, Ef, vf), some are specific to a particular type of the
+ PulloutModel.
+ """
+
+
+class CBMConstantBond(CrackBridgeModel):
+ """
+ Return the matrix stress profile of a crack bridge for a given control slip
+ at the loaded end
+ """
+ Em = bu.Float(25e3, MAT=True) # [MPa] matrix modulus
+ Ef = bu.Float(180e3, MAT=True) # [MPa] fiber modulus
+ vf = bu.Float(0.01, CS=True) # [-] reinforcement ratio
+ T = bu.Float(12., MAT=True) # [N/mm^3] bond intensity
+
+ ipw_view = bu.View(
+ bu.Item('Em'),
+ bu.Item('Ef'),
+ bu.Item('vf'),
+ bu.Item('T'),
+ )
+
+class CrackBridgeRespSurface(bu.Model):
+
+ cb = bu.Instance(CrackBridgeModel)
+
+ ## Crack bridge with constant bond
+ def get_sig_m(self, z, sig_c): # matrix stress (*\label{sig_m}*)
+ T = self.cb.T
+ vf = self.cb.vf
+ Em = self.cb.Em
+ Ef = self.cb.Ef
+ sig_m = np.minimum(z * T * vf / (1 - vf), Em * sig_c / (vf * Ef + (1 - vf) * Em))
+ return sig_m
+
+ def get_eps_f(self, z, sig_c): # reinforcement strain (*\label{sig_f}*)
+ vf = self.cb.vf
+ Ef = self.cb.Ef
+ sig_m = self.get_sig_m(z, sig_c)
+ eps_f = (sig_c - sig_m * (1 - vf)) / vf / Ef
+ return eps_f
+
+ sig_c_slider: float = bu.Float(1.0, BC=True)
+
+ ipw_view = bu.View(
+ bu.Item('sig_c_slider')
+ )
+
+ @staticmethod
+ def subplots(fig):
+ ax = fig.subplots(1,1)
+ ax1 = ax.twinx()
+ return ax, ax1
+
+ def update_plot(self, axes):
+ ax, ax1 = axes
+ x_range = np.linspace(-300,300,1000)
+ z_range = np.abs(x_range)
+
+ sig_m_range = self.get_sig_m(z_range, self.sig_c_slider)
+ eps_f_range = self.get_eps_f(z_range, self.sig_c_slider)
+ sig_max = np.max(sig_m_range)
+ eps_max = np.max(eps_f_range)
+
+ ax.plot(x_range, sig_m_range, color='black')
+ ax.fill_between(x_range, sig_m_range, color='gray', alpha=0.1)
+ ax.set_ylim(ymin=-0.03*sig_max)
+ ax1.plot(x_range, eps_f_range, color='blue')
+ ax1.fill_between(x_range, eps_f_range, color='blue', alpha=0.1)
+ ax1.set_ylim(ymin=-0.03*eps_max)
+
+
+class PMCM(bu.Model):
+ name='Fragmentation'
+ sig_cu = bu.Float(10.0, MAT=True) # [MPa] composite strength
+ sig_mu = bu.Float(3.0, MAT=True) # [MPa] matrix strength
+ m = bu.Float(4, MAT=True) # Weibull shape modulus
+ ## Crack tracing algorithm
+ n_x = bu.Int(200, MAT=True)
+ L_x = bu.Float(380, MAT=True)
+
+ cb = bu.Instance(CrackBridgeModel)
+ def _cb_default(self):
+ return CBMConstantBond()
+
+ cb_rs = tr.Property
+ @tr.cached_property
+ def _get_cb_rs(self):
+ return CrackBridgeRespSurface(cb=self.cb)
+ #
+ # # %%
+ #
+ # cb_rs = bu.Instance(CrackBridgeRespSurface, ())
+
+ tree = ['cb', 'cb_rs']
+
+ ipw_view = bu.View(
+ bu.Item('n_x'),
+ bu.Item('L_x'),
+ bu.Item('sig_cu'),
+ bu.Item('sig_mu'),
+ bu.Item('m'),
+ )
+
+
+ ## Specimen discretization
+ def get_z_x(self, x, XK): # distance to the closest crack (*\label{get_z_x}*)
+ z_grid = np.abs(x[:, np.newaxis] - np.array(XK)[np.newaxis, :])
+ return np.amin(z_grid, axis=1)
+
+ def get_sig_c_z(self, sig_mu, z, sig_c_pre):
+ # crack initiating load at a material element
+ fun = lambda sig_c: sig_mu - self.cb_rs.get_sig_m(z, sig_c)
+ try: # search for the local crack load level
+ return newton(fun, sig_c_pre)
+ except (RuntimeWarning, RuntimeError):
+ # solution not found (shielded zone) return the ultimate composite strength
+ return self.sig_cu
+
+ def get_sig_c_K(self, z_x, x, sig_c_pre, sig_mu_x):
+ # crack initiating loads over the whole specimen
+ get_sig_c_x = np.vectorize(self.get_sig_c_z)
+ sig_c_x = get_sig_c_x(sig_mu_x, z_x, sig_c_pre)
+ y_idx = np.argmin(sig_c_x)
+ return sig_c_x[y_idx], x[y_idx]
+
+ cracking_history = tr.Property(depends_on='state_changed')
+ @tr.cached_property
+ def _get_cracking_history(self):
+ return self.get_cracking_history()
+
+ def get_cracking_history(self, update_progress=None):
+ L_x, n_x = self.L_x, self.n_x
+ sig_mu = self.sig_mu
+ m = self.m
+ x = np.linspace(0, L_x, n_x) # specimen discretization (*\label{discrete}*)
+ sig_mu_x = sig_mu * np.random.weibull(m, size=n_x) # matrix strength (*\label{m_strength}*)
+
+ vf = self.cb_rs.cb.vf
+ Em = self.cb_rs.cb.Em
+ Ef = self.cb_rs.cb.Ef
+
+ Ec = Em * (1 - vf) + Ef * vf # [MPa] mixture rule
+
+ XK = [] # recording the crack postions
+ sig_c_K = [0.] # recording the crack initating loads
+ eps_c_K = [0.] # recording the composite strains
+ CS = [L_x, L_x / 2] # crack spacing
+ sig_m_x_K = [np.zeros_like(x)] # stress profiles for crack states
+
+ idx_0 = np.argmin(sig_mu_x)
+ XK.append(x[idx_0]) # position of the first crack
+ sig_c_0 = sig_mu_x[idx_0] * Ec / Em
+ sig_c_K.append(sig_c_0)
+ eps_c_K.append(sig_mu_x[idx_0] / Em)
+
+ while True:
+ z_x = self.get_z_x(x, XK) # distances to the nearest crack
+ sig_m_x_K.append(self.cb_rs.get_sig_m(z_x, sig_c_K[-1])) # matrix stress
+ sig_c_k, y_i = self.get_sig_c_K(z_x, x, sig_c_K[-1], sig_mu_x) # identify next crack
+ if sig_c_k == self.sig_cu: # (*\label{no_crack}*)
+ break
+ if update_progress: # callback to user interface
+ update_progress(sig_c_k)
+ XK.append(y_i) # record crack position
+ sig_c_K.append(sig_c_k) # corresponding composite stress
+ eps_c_K.append( # composite strain - integrate the strain field
+ np.trapz(self.cb_rs.get_eps_f(self.get_z_x(x, XK), sig_c_k), x) / np.amax(x)) # (*\label{imple_avg_strain}*)
+ XK_arr = np.hstack([[0], np.sort(np.array(XK)), [L_x]])
+ CS.append(np.average(XK_arr[1:] - XK_arr[:-1])) # crack spacing
+
+ sig_c_K.append(self.sig_cu) # the ultimate state
+ eps_c_K.append(np.trapz(self.cb_rs.get_eps_f(self.get_z_x(x, XK), self.sig_cu), x) / np.amax(x))
+ CS.append(CS[-1])
+ if update_progress:
+ update_progress(sig_c_k)
+ return np.array(sig_c_K), np.array(eps_c_K), sig_mu_x, x, np.array(CS), np.array(sig_m_x_K)
+
+ @staticmethod
+ def subplots(fig):
+ ax1, ax2 = fig.subplots(1,2)
+ ax11 = ax1.twinx()
+ return ax1, ax11, ax2
+
+ def plot(self, axes):
+ ax, ax_cs, ax_sig_x = axes
+ sig_c_K, eps_c_K, sig_mu_x, x, CS, sig_m_x_K = self.cracking_history
+ #sig_c_K, eps_c_K, sig_mu_x, x, CS, sig_m_x_K, sig_m_x_K1 = self.cracking_history
+ n_c = len(eps_c_K) - 2 # numer of cracks
+ ax.plot(eps_c_K, sig_c_K, marker='o', label='%d cracks:' % n_c)
+ ax.set_xlabel(r'$\varepsilon_\mathrm{c}$ [-]');
+ ax.set_ylabel(r'$\sigma_\mathrm{c}$ [MPa]')
+ ax_sig_x.plot(x, sig_mu_x, color='orange')
+ ax_sig_x.fill_between(x, sig_mu_x, 0, color='orange', alpha=0.1)
+ ax_sig_x.set_xlabel(r'$x$ [mm]');
+ ax_sig_x.set_ylabel(r'$\sigma$ [MPa]')
+ ax.legend()
+ eps_c_KK = np.array([eps_c_K[:-1], eps_c_K[1:]]).T.flatten()
+ CS_KK = np.array([CS[:-1], CS[:-1]]).T.flatten()
+ ax_cs.plot(eps_c_KK, CS_KK, color='gray')
+ ax_cs.fill_between(eps_c_KK, CS_KK, color='gray', alpha=0.2)
+ ax_cs.set_ylabel(r'$\ell_\mathrm{cs}$ [mm]');
+
+ def update_plot(self, axes):
+ self.plot(axes)
+
+# %%
+
+def render_widgets():
+ ## Interactive application
+ import ipywidgets as ipw
+
+ n_steps = 20
+ margs_sliders = {
+ name: ipw.FloatSlider(description=desc, value=val,
+ min=minval, max=maxval, step=(maxval - minval) / n_steps,
+ continuous_update=False)
+ for name, desc, val, minval, maxval in [
+ ('Em', r'\(E_\mathrm{m}\;\)[MPa]', 28e3, 1e3, 50e3),
+ ('Ef', r'\(E_\mathrm{f}\;\)[MPa]', 180e3, 1e3, 250e3),
+ ('vf', r'\(V_\mathrm{f}\;\)[-]', 0.01, 0.00001, 0.4),
+ ('T', r'\(T\;\)[N/mm$^3$]', 8, 0.0001, 20),
+ ('sig_cu', r'\(\sigma_\mathrm{cu}\;\)[MPa]', 10, 3, 100),
+ ('sig_mu', r'\(\sigma_\mathrm{mu}\;\)[MPa]', 5.0, 1, 10),
+ ('m', r'\(m\;\)[-]', 4, 0.8, 100),
+ ('L_x', r'\(L\;\)[mm]', 500, 200, 2000)
+ ]
+ }
+ pmcm_args = ['sig_cu', 'sig_mu', 'm', 'L_x', 'n_x']
+ cb_args = ['Em', 'Ef', 'vf', 'T']
+
+ margs_sliders['n_x'] = ipw.IntSlider(description='$n_\mathrm{points}$', value=200,
+ min=20, max=1000, step=10, continuous_update=False)
+ crack_slider = ipw.IntSlider(description='crack', value=0, min=0, max=1, step=1)
+ progress = ipw.FloatProgress(min=0, max=1) # instantiate the bar
+
+ fig, (ax, ax_sig_x) = plt.subplots(1, 2, figsize=(8, 3), tight_layout=True)
+ ax_cs = ax.twinx()
+
+
+ def update_progress(sig):
+ progress.value = sig
+
+
+ pmcm = PMCM()
+
+ def init():
+ for key, sl in margs_sliders.items():
+ globals()[key] = sl.value
+ sig_c_K, eps_c_K, sig_mu_x, x, CS, sig_m_x_K = pmcm.get_cracking_history(update_progress) # (*\label{calc_curve}*)
+ progress.max = margs_sliders['sig_cu'].value
+ ax.plot(eps_c_K, sig_c_K, marker='o') # (*\label{show_curve1}*)
+ ax_sig_x.plot(x, sig_mu_x, color='red')
+
+
+ current_sig_m_x_K = []
+ current_x = []
+
+
+ def reset_crack_slider(x, sig_m_x_K):
+ global current_sig_m_x_K, current_x, sig_m_line, sig_eps_marker
+ current_sig_m_x_K = sig_m_x_K
+ current_x = x
+ n_cracks = len(sig_m_x_K)
+ crack_slider.max = n_cracks - 1
+ crack_slider.value = 0
+ sig_m_line, = ax_sig_x.plot(x, sig_m_x_K[0])
+ sig_eps_marker, = ax.plot([0], [0], color='magenta', marker='o')
+
+
+ def update_crack_slider(crack):
+ global sig_m_line, sig_eps_marker
+ global sig_c_K, eps_c_K
+ if len(current_sig_m_x_K) > 0:
+ sig_m_line.set_ydata(current_sig_m_x_K[crack])
+ sig_eps_marker.set_data(eps_c_K[crack], sig_c_K[crack])
+
+
+ def update(**mparams):
+ global sig_c_K, eps_c_K
+ for key, val in mparams.items():
+ globals()[key] = val
+ ax.clear()
+ ax_cs.clear()
+ ax_sig_x.clear()
+ ax.set_title('stress-strain and crack spacing')
+ ax_sig_x.set_title('matrix strength and stress profiles')
+ sig_c_K, eps_c_K, sig_mu_x, x, CS, sig_m_x_K = pmcm.get_cracking_history(update_progress) # (*\label{calc_curve}*)
+ n_c = len(eps_c_K) - 2 # numer of cracks
+ ax.plot(eps_c_K, sig_c_K, marker='o', label='%d cracks:' % n_c) # (*\label{show_curve1}*)
+ ax.set_xlabel(r'composite strain $\varepsilon_\mathrm{c}$ [-]');
+ ax.set_ylabel(r'composite stress $\sigma_\mathrm{c}$ [MPa]')
+ ax_sig_x.plot(x, sig_mu_x, color='orange')
+ ax_sig_x.fill_between(x, sig_mu_x, 0, color='orange', alpha=0.1)
+ ax_sig_x.set_xlabel(r'$x$ [mm]');
+ ax_sig_x.set_ylabel(r'matrix stress $\sigma_\mathrm{m}$ [MPa]')
+ ax.legend()
+ eps_c_KK = np.array([eps_c_K[:-1], eps_c_K[1:]]).T.flatten()
+ CS_KK = np.array([CS[:-1], CS[:-1]]).T.flatten()
+ ax_cs.plot(eps_c_KK, CS_KK, color='gray')
+ ax_cs.fill_between(eps_c_KK, CS_KK, color='gray', alpha=0.2)
+ ax_cs.set_ylabel(r'crack spacing $\ell_\mathrm{cs}$ [mm]');
+ reset_crack_slider(x, sig_m_x_K)
+
+
+ def slider_layout():
+ hbox_pr = ipw.HBox([progress])
+ hbox = ipw.HBox([crack_slider])
+ pmcm_layout = ipw.Layout(grid_template_columns='1fr')
+ pmcm_list = [ipw.Label('Crack tracing parameters')] + \
+ [margs_sliders[arg] for arg in pmcm_args]
+ pmcm_col = ipw.GridBox(pmcm_list, layout=pmcm_layout)
+ cb_layout = ipw.Layout(grid_template_columns='1fr')
+ cb_list = [ipw.Label('Crack bridge parameters')] + \
+ [margs_sliders[arg] for arg in cb_args]
+ cb_col = ipw.GridBox(cb_list, layout=cb_layout)
+ alg_params = ipw.HBox([pmcm_col, cb_col])
+ box = ipw.VBox([hbox_pr, hbox, alg_params])
+ display(box)
+
+
+ init()
+ slider_layout()
+ ipw.interactive_output(update_crack_slider, {'crack': crack_slider})
+ ipw.interactive_output(update, margs_sliders);