<jupyter_start><jupyter_text>Define the flow potential<jupyter_code>import sympy as sp
sp.init_printing()<jupyter_output><empty_output><jupyter_text>Logic behind the flow potential extensionTo account for a broad range of applications we proposean extension of the flow potential with an adjustable level of interaction between the damage in the normaland tangential directions. Recall that the evolution equationsare defined as\begin{align} \dot{\omega}_\mathrm{N} &= \lambda \frac{\partial \varphi}{ \partial {\omega}_\mathrm{N} } \\ \dot{\omega}_\mathrm{T} &= \lambda \frac{\partial \varphi}{ \partial {\omega}_\mathrm{T} }\end{align}If the damage due to normal direction directly affects the tangential version we assumethat also the opposite affect is the same. The respective derivatives of the potential must be equal, i.e.\begin{align}\frac{\partial \varphi}{ \partial {\omega}_\mathrm{N} } = \frac{\partial \varphi}{ \partial {\omega}_\mathrm{T} }\end{align}<jupyter_code>S_N, S_T, eta = sp.symbols('S_N, S_T, eta')
Y_N, Y_T = sp.symbols('Y_N, Y_T')
omega_N, omega_T = sp.symbols('omega_N, omega_T')
c_N, c_T = sp.symbols('c_N, c_T')
phi_N_ext = (  # (1-omega_N)**c_N * (
    (1*Y_N*(Y_N + eta * Y_T)) /
    (2*(S_N - eta * (S_N - (S_N + S_T)/2)))
)
phi_T_ext = (  # (1-omega_T)**c_T * (
    (1*Y_T*(Y_T + eta * Y_N)) /
    (2*(S_T - eta * (S_T - (S_N + S_T)/2)))
)
phi_ext = sp.simplify(phi_N_ext) + sp.simplify(phi_T_ext)
phi_ext
phi_ext.subs(eta,1)
sp.simplify(phi_ext.diff(Y_N).subs(eta,1))
sp.simplify(phi_ext.diff(Y_T).subs(eta,1))
sp.simplify(phi_ext.diff(Y_N).subs(eta,0))<jupyter_output><empty_output><jupyter_text>Backward derivation of the potential from desired flow direction Find a potential function \begin{align}\varphi = a Y_\mathrm{N}^2 + \eta b Y_\mathrm{N}(Y_\mathrm{N} + Y_\mathrm{T})+ c Y_\mathrm{T}^2 + \eta d Y_\mathrm{T}(Y_\mathrm{T} + Y_\mathrm{N})\end{align}which fulfills\begin{align}\left. \frac{\partial \varphi}{\partial Y_\mathrm{N}}\right|_{\eta=0} &= (1 - \omega_\mathrm{N})^{c_\mathrm{N}} \frac{Y_\mathrm{N}}{S_\mathrm{N}} \\\left. \frac{\partial \varphi}{\partial Y_\mathrm{T}}\right|_{\eta=0} &= (1 - \omega_\mathrm{T})^{c_\mathrm{T}} \frac{Y_\mathrm{T}}{S_\mathrm{T}} \\\left. \frac{\partial \varphi}{\partial Y_\mathrm{N}}\right|_{\eta=1} = \left. \frac{\partial \varphi}{\partial Y_\mathrm{T}}\right|_{\eta=1} &= \left(1 - \frac{1}{2}(\omega_\mathrm{N} + \omega_\mathrm{T})\right)^{\frac{c_\mathrm{N}+c_\mathrm{T}}{2}} \frac{Y_\mathrm{N} + Y_\mathrm{T}}{S_\mathrm{N} + S_\mathrm{T}}\end{align}<jupyter_code>a, b, c, d = sp.symbols('a,b,c,d')
H = lambda x: sp.Piecewise( (0, x <=0 ), (1, True) )
H_switch, sig_pi = sp.symbols(r'H(\sigma^{\pi}), \sigma^\pi', real=True)
phi_ext = a * Y_N**2 + b * eta * Y_N*(Y_N + Y_T) + c * Y_T**2 + d * eta * Y_T*(Y_T+Y_N)
d_phi_N_0 = phi_ext.diff(Y_N).subs(eta,0)
a_solved = sp.solve(sp.Eq(d_phi_N_0, (1 - omega_N)**c_N * Y_N / S_N * H_switch ), a )[0]
#a_solved = sp.solve( sp.Eq( d_phi_N_0, (1 - omega_N)**c_N * Y_N / S_N ), a )[0]
d_phi_N_0, a_solved
d_phi_T_0 = phi_ext.diff(Y_T).subs(eta,0)
c_solved = sp.solve( sp.Eq( d_phi_T_0, (1 - omega_T)**c_T * Y_T / S_T ), c )[0]
d_phi_T_0, c_solved
phi_ext_ac = phi_ext.subs({a: a_solved, c: c_solved})
phi_ext_ac
d_phi_N_1 = phi_ext_ac.diff(Y_N).subs(eta,1)
d_phi_N_1
d_phi_T_1 = phi_ext_ac.diff(Y_T).subs(eta,1)
d_phi_T_1
d_phi_1_req = (1 - (omega_N + omega_T)/2)**((c_N+c_T)/2) * (Y_N + Y_T) / (S_N + S_T) 
d_phi_1_req
bd_solved = sp.solve({sp.Eq(d_phi_N_1, d_phi_1_req), sp.Eq(d_phi_T_1, d_phi_1_req)},[b,d])
sp.simplify(bd_solved[b])
phi_abcd = phi_ext_ac.subs(bd_solved)
sp.simplify(phi_abcd)
sp.simplify(phi_abcd.diff(Y_T).subs(eta,1))
sp.simplify(phi_abcd.diff(Y_N).subs(eta,sp.Rational(1,2)))<jupyter_output><empty_output><jupyter_text>Geometric mean<jupyter_code>a,b,c,d = sp.symbols('a,b,c,d')
H = lambda x: sp.Piecewise( (0, x <=0 ), (1, True) )
H_switch, sig_pi = sp.symbols(r'H(\sigma^{\pi}), \sigma^\pi', real=True)
phi2_ext = a * Y_N**2 + b * eta * Y_N*(Y_N+Y_T) + c * Y_T**2 + d * eta * Y_T*(Y_N+Y_T)
phi2_ext
d_phi2_N_0 = phi2_ext.diff(Y_N).subs(eta,0)
a2_solved = sp.solve( sp.Eq( d_phi2_N_0, (1 - omega_N)**c_N * Y_N / S_N * H_switch ), a )[0]
#a_solved = sp.solve( sp.Eq( d_phi_N_0, (1 - omega_N)**c_N * Y_N / S_N ), a )[0]
d_phi2_N_0, a2_solved
d_phi2_T_0 = phi2_ext.diff(Y_T).subs(eta,0)
c2_solved = sp.solve( sp.Eq( d_phi2_T_0, (1 - omega_T)**c_T * Y_T / S_T ), c )[0]
d_phi2_T_0, c2_solved
phi2_ext_ac = phi2_ext.subs({a: a2_solved, c: c2_solved})
phi2_ext_ac
d_phi2_N_1 = phi2_ext_ac.diff(Y_N).subs(eta,1)
d_phi2_N_1
d_phi2_T_1 = phi2_ext_ac.diff(Y_T).subs(eta,1)
d_phi2_T_1
c_NT = sp.sqrt(c_N*c_T)
S_NT = sp.sqrt(S_N*S_T)
d_phi_2_req = (1 - sp.sqrt(omega_N*omega_T))**(c_NT) * (Y_N+Y_T) / (2*S_NT) 
d_phi_2_req
bd2_solved = sp.solve({sp.Eq(d_phi2_N_1, d_phi_2_req), sp.Eq(d_phi2_T_1, d_phi_2_req)},[b,d])
sp.simplify(bd2_solved[b])
phi2_abcd = phi2_ext_ac.subs(bd2_solved)
sp.simplify(phi2_abcd).subs(eta,0)
omega_T_solved_0 = sp.simplify(phi2_abcd.diff(Y_T).subs(eta,0))
omega_T_solved_0
omega_N_solved_0 = sp.simplify(phi2_abcd.diff(Y_N).subs(eta,0))
omega_N_solved_0
omega_N_solved_1 = sp.simplify(phi2_abcd.diff(Y_N).subs(eta,1))
omega_N_solved_1
omega_T_solved_1 = sp.simplify(phi2_abcd.diff(Y_T).subs(eta,1))
omega_T_solved_1<jupyter_output><empty_output><jupyter_text>**Archived compact version:** derived by Mirek for the case of dropped fundamental damage terms $(1 - \omega_\mathrm{N})^{c_\mathrm{N}}$ and $(1 - \omega_\mathrm{T})^{c_\mathrm{T}}$<jupyter_code>phi_simplified = (
    (( 1 - eta ) * ( Y_T**2 * S_N**2 + Y_N**2 * S_T**2 ) +
    S_N * S_T * ( Y_T**2 + Y_N**2 + 2 * Y_T * Y_N * eta )) /
    (2* S_N * S_T * ( S_N + S_T ))
)
phi_simplified<jupyter_output><empty_output><jupyter_text>Scaling of phi between coupled and uncoupled cases<jupyter_code>def avg(var1, var2):
    return (var1 + var2) / 2
r = sp.symbols(r'r', positive=True)
c_NT = avg(c_N, c_T)
S_NT = avg(S_N, S_T)
omega_NT = avg(omega_N, omega_T)
Y_NT = avg(Y_N, Y_T)
phi_N = (1 - omega_N)**(c_N) * S_N/(r+1) * (Y_N/S_N)**(r+1) # * H_switch
phi_T = (1 - omega_T)**(c_T) * S_T/(r+1) * (Y_T/S_T)**(r+1)
phi_NT  = (1 - omega_NT)**(c_NT) * 2*S_NT/(r+1) * (Y_NT/S_NT)**(r+1)
phi_ = (1 - eta) * (phi_N + phi_T) + eta * phi_NT
sp.simplify(phi_).subs(r,1).subs(eta,0).diff(Y_T)
sp.simplify(phi_).subs(r,1).subs(eta,0).diff(Y_N)<jupyter_output><empty_output><jupyter_text>Derivation of potential from a damage function<jupyter_code>A_d, Y_N, R_N, lambda_N = sp.symbols('A_d, Y_N, R_N, lambda_N')
dot_Y_N, dot_R_N = sp.symbols(r'\dot{Y}_N, \dot{R}_N')
E_N , eps_N = sp.symbols(r'E_N, \epsilon_N')
dot_eps_N = sp.symbols(r'\dot{\epsilon}_N')
r_N, dot_r_N = sp.symbols(r'r_N, \dot{r}_N')
f_N = sp.Function('f_N')(r_N)
omega_N = sp.symbols(r'\omega_N')
dot_omega_N = sp.symbols(r'\dot{\omega}_N')
psi_omega = sp.Rational(1,2)*(1 - omega_N)*E_N*eps_N**2 + f_N
psi_omega
sigma_N_ = psi_omega.diff(eps_N)
Y_N_ = -psi_omega.diff(omega_N)
R_N_ = psi_omega.diff(r_N)
Y_N_
dot_R_N_ = R_N_.diff(r_N) * dot_r_N
dot_R_N_
dot_Y_N_ = Y_N_.diff(eps_N) * dot_eps_N + Y_N_.diff(omega_N) * dot_omega_N
dot_Y_N_
f_omega = Y_N - R_N
f_omega
dot_r_N_ = lambda_N * f_omega.diff(R_N)
dot_r_N_
dot_omega_N_ = lambda_N * f_omega.diff(Y_N)
dot_omega_N_
dot_f_omega_ = f_omega.diff(Y_N) * dot_Y_N_ + f_omega.diff(R_N) * dot_R_N_
dot_f_omega_
dot_f_omega_lambda_ = dot_f_omega_.subs(dot_r_N, dot_r_N_)
lambda_N_solved = sp.solve(dot_f_omega_lambda_, lambda_N)[0]
lambda_N_solved
dot_omega_N_solved = dot_omega_N_.subs(lambda_N, lambda_N_solved)
dot_omega_N_solved<jupyter_output><empty_output><jupyter_text>Damage function from MS1 - N **reference [74]**<jupyter_code>omega_N_paper = 1 - 1 / (1 + A_d * Y_N_) 
omega_N_paper
dot_omega_N_paper = omega_N_paper.diff(eps_N) * dot_eps_N
dot_omega_N_paper<jupyter_output><empty_output><jupyter_text>**from potential**<jupyter_code>f_prime_r_N = 1 / A_d * -r_N / (1 + r_N)
f_prime_r_N
f_prime_r_N.diff(r_N)
dot_omega_N_potential = sp.simplify(- E_N * eps_N * dot_eps_N / f_prime_r_N.diff(r_N))
dot_omega_N_potential
r_N_solved = sp.solve( sp.Eq( dot_omega_N_paper, dot_omega_N_potential ), r_N )
r_N_solved
dot_omega_N_potential_solved = sp.simplify(dot_omega_N_potential.subs(r_N, r_N_solved[1]))
dot_omega_N_potential_solved
sp.simplify(dot_omega_N_paper)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18343337334933973, 'fraction_numerical': 0.022929171668667467, 'mean_word_length': 5.648842777334398, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '521330', 'license': 'MIT', 'n_tokens_mistral': 4172, 'n_tokens_neox': 4091, 'n_words': 624, 'path': 'notebooks/flow_potential.ipynb', 'repo_name': 'bmcs-group/bmcs_slide', 'size': '8330'}