You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
# Here the ordinary differential equations (ODEs), defined by $f$ describes the electro-chemical state of the muscle cells, characterized by the state vector $s$. Furthermore v is the transmembrane potential, $\lambda$ is the fiber stretch ratio and $\dot{\lambda} = \frac{\partial \lambda}{\partial t}$ is the stretch rate. The electrical state of the tissue is described by the Bi-domain model, which could also be reduced to the Mondomain model (see {cite}`sundnes2007computing` for more info). The function $I_{\text{ion}}$ describes the ionic current accross the cell membrane, while $M_i$ and $M_e$ are intracellular and extracellular tissue conductivities respectively.
29
+
# Here the ordinary differential equations (ODEs), defined by $f$ describes the electro-chemical state of the muscle cells, characterized by the state vector $s$. Furthermore $v$ is the transmembrane potential, $\lambda$ is the fiber stretch ratio and $\dot{\lambda} = \frac{\partial \lambda}{\partial t}$ is the stretch rate. The electrical state of the tissue is described by the Bi-domain model, which could also be reduced to the Mondomain model (see {cite}`sundnes2007computing` for more info). The function $I_{\text{ion}}$ describes the ionic current accross the cell membrane, while $M_i$ and $M_e$ are intracellular and extracellular tissue conductivities respectively.
30
30
#
31
31
# The mechanical equilibrium are described by the following system of equations
# Here $\Psi$ is a strain energy function which relates stress to strain. For examples we can choose a NeoHookean material model (which is typically for rubberlike material)
52
+
# Here $\Psi$ is a strain energy function which relates stress to strain and depends on the material under consideration (This is often denoted as the "constitutive relation"). For examples we can choose a NeoHookean material model (which is typical for rubber-like materials)
# The active stress $\mathbf{P}_a$ depends on the electro-chemical state, and is typically described via a model for the crossbrigde dynamics.
59
59
#
60
-
# Now the mechanics system is typically solved for the displacement field $\mathbf{u}$, which relates coordinates in the reference configuration ($\mathbf{X}$) to coordinates in the current configuration ($\mathbf{x}$), by $\mathbf{u} = \mathbf{x} - \mathbf{X}$. We also have the deformation gradient $\mathbf{F} = \nabla \mathbf{u} + \mathbf{I}$. Furthermore the stretch is typically defined as the stretch along the fiber direction. If we denote the fiber orientation in the reference configuration as $\mathbf{f}_0$, then $\lambda = | \mathbf{F}\mathbf{f}_0 | = | \mathbf{f} | $. For example one can consider a tissue slab with fibers oriented in the $x$-direction. In this case $\mathbf{f}_0 = (1 \; 0 \; 0)^T$
60
+
# Now the mechanics system is typically solved for the displacement field $\mathbf{u}$, which relates coordinates in the reference configuration ($\mathbf{X}$) to coordinates in the current configuration ($\mathbf{x}$), by $\mathbf{u} = \mathbf{x} - \mathbf{X}$. We also have the deformation gradient $\mathbf{F} = \nabla \mathbf{u} + \mathbf{I}$. Furthermore the stretch is typically defined as the stretch along the fiber direction. If we denote the fiber orientation in the reference configuration as $\mathbf{f}_0$, then $\lambda = | \mathbf{F}\mathbf{f}_0 | = | \mathbf{f} | $. For example one can consider a tissue slab with fibers oriented in the $x$-direction, in which case we would have $\mathbf{f}_0 = (1 \; 0 \; 0)^T$
61
61
#
62
62
# ## Reduction to 0D
63
63
#
64
-
# To reduce the problem to 0D, one need to consider a simplified experiment. For example if we consider a rectangle length $l_x$, heigth $l_y$ and depth $l_z$. Now imaging that we apply a deformation that changes the length from $l_x$ to $l_x + \Delta x = L_x$. We can quantify this deformation using the stretch
64
+
# To reduce the problem to 0D, one need to consider a simplified experiment. For example we can consider a rectangle of length $l_x$, a heigth of $l_y$ and a depth of $l_z$. Now imaging that we apply a deformation that changes the length from $l_x$ to $l_x + \Delta x = L_x$. We can quantify this deformation using the stretch
# If we now assume the material is incompressible, and that the compression along the $y$- and $z$-axis are equal then we can write the following deformation gradient
70
+
# If we now assume that the material is incompressible, and that the compression along the $y$- and $z$-axis are equal then we can write the following deformation gradient
# Note here that the variable $p$ is introduced. This variable acts as a Lagrange multiplier to enforce incompressibilty
155
+
# Note here that the variable $p$ is introduced. This variable acts as a Lagrange multiplier to enforce incompressibility
156
156
157
157
print(mech_model.compressibility)
158
158
print(mech_model.compressibility.str())
@@ -165,7 +165,7 @@
165
165
# \mathbf{P}_a = \mathbf{P}_a(s)
166
166
# $$
167
167
#
168
-
# We will now update our model to use an active stress model, which contains a variable $T_a$ that describes the stregth of the active force and a fiber direction (which in our case is chosen to be along the $x$-axis)
168
+
# We will now update our model to use an active stress model, which contains a variable $T_a$ that describes the strength of the active force and a fiber direction (which in our case is chosen to be along the $x$-axis)
169
169
170
170
active_model=zero_mech.active.ActiveStress()
171
171
display(active_model)
@@ -223,7 +223,7 @@
223
223
mon=model["monitor_values"]
224
224
# -
225
225
226
-
# We set the initial $\lambda$ to 1.0 and $\dot{\lambda}$ to 0.0. Since we would need to update $\lambda$ and $\dot{\lambda}$ we need to also keep track of the previous value of $\lambda$.
226
+
# We set the initial $\lambda$ to 1.0 and $\dot{\lambda}$ to 0.0. Since we would need to update $\lambda$ and $\dot{\lambda}$, we need to also keep track of the previous value of $\lambda$.
227
227
228
228
# +
229
229
lmbda_index=model["parameter_index"]("lmbda")
@@ -237,7 +237,7 @@
237
237
238
238
# -
239
239
240
-
# Now, we will create a function that will solve the mechanics problem (i.e P_xx = P_yy = 0) for $\lambda$ and $p$, given an active stress $T_a$
240
+
# Now, we will create a function that will solve the mechanics problem (i.e $P_{xx} = P_{yy} = 0$) for $\lambda$ and $p$, given an active stress $T_a$
241
241
242
242
deffunc(x, Ta):
243
243
lmbda, p=x
@@ -254,7 +254,7 @@ def func(x, Ta):
254
254
returnnp.array([P11, P22], dtype=np.float64)
255
255
256
256
257
-
# Now let us loop over the time steps and in each iteration we first solve the EP problem by calling the `fgr` (forward generalized rush larsen) function, and then we solve the mechanics problem using a root-finding algorithm. An important observation to make here is that we do not update $\lambda$ or $\dot{\lambda}$ here (we have commented out the relevant lines that would update these values in the EP model)
257
+
# Now let us loop over the time steps and in each iteration we first solve the EP problem by calling the `fgr` (forward generalized rush larsen) function, and then we solve the mechanics problem using a root-finding algorithm. An important observation to make here is that we do not update $\lambda$ or $\dot{\lambda}$ (we have commented out the relevant lines that would update these values in the EP model)
258
258
259
259
# +
260
260
fromscipy.optimizeimportroot
@@ -385,9 +385,9 @@ def func(x, Ta):
385
385
plt.show()
386
386
387
387
388
-
# Yikes!! That did not work well. We see that the solution starts oscillating. The reason this happens is because $T_a$ depends heavily on $\lambda$ and $\dot{\lambda}$, but we also have that $\lambda$ and $\dot{\lambda}$ depends heavily on $T_a$. For example if we stretch the material (i.e increase $\lambda$), then the [Frank Starling law](https://en.wikipedia.org/wiki/Frank%E2%80%93Starling_law) states that $T_a$ should be increased, but when $T_a$ increases we get more force that compress the material which will bring $\lambda$ down and it is not hard do see why this system can start oscillating.
388
+
# Yikes!! That did not work well. We see that the solution starts oscillating. The reason this happens is because $T_a$ depends heavily on $\lambda$ and $\dot{\lambda}$, but we also have that $\lambda$ and $\dot{\lambda}$ depends heavily on $T_a$. For example if we stretch the material (i.e increase $\lambda$), then the [Frank Starling law](https://en.wikipedia.org/wiki/Frank%E2%80%93Starling_law) states that $T_a$ should be increased, but when $T_a$ increases we get more force that compress the material which will bring $\lambda$ down, and it is not hard do see why this system can start oscillating.
389
389
#
390
-
# The solution to the problem is to make sure we update $T_a$ whenever we solve the mechanics model, which means that we would need to run the ODE systen during every iteration whithin the root finding algorithm. In other words, let us update the function as follows
390
+
# The solution to the problem is to make sure we update $T_a$ whenever we solve the mechanics model, which means that we would need to integrate the ODE systen during every iteration within the root finding algorithm. In other words, let us update the function as follows
# Note that whenever we now get a new guess for $\lambda$ (and $p$), we update $T_a$ accordingly. Also now that we also need to pass in the previous $\lambda$ (in order to compute $\dot{\lambda}$ as well as the previous state, the parameters and the time step. Our time loop will therefore need to be updated accordingly.
415
+
# Note that whenever we now get a new guess for $\lambda$ (and $p$), we update $T_a$ accordingly. We now need to also pass the previous $\lambda$ (in order to compute $\dot{\lambda}$) as well as the previous state, the parameters and the time step. Our time loop will therefore need to be updated accordingly.
0 commit comments