notebook.community

Edit and run

This notebook contains all codes necessary to reproduce the study published in the article "A Characteristic-based constitutive law for dispersed fibers". The notebook was created using SageMathCloud (SMC) (http://cloud.sagemath.com), which is the best place to run and modify the current code.

To run the code on SMC, simply copy this notebook over to SMC. Open the notebook and set the jupyter notebook kernel to SageMath. All source codes were written in Python and the computational engine is SageMath. Here is a tutorial for SageMath: http://doc.sagemath.org/html/en/tutorial/index.html

Set the kernel to SageMath

Execute all codes in this notebook

Copyright © Liang Ge (liang.ge@gmail.com)

In this analysis we assumed the mean fiber direction is along the $x$ axis. Two Eulerian angles $\Theta$ and $\Phi$ define the fiber direction vector, $e_f=(\cos(\Theta), \sin(\Theta)\cos(\Phi), \sin(\Theta)\sin(\Phi))$.

Fig. 1 Fiber direction vector $e_f$ defined by two Eulerian angles $\Theta$ and $\Phi$

The fiber stretching ratio after deformation is $\lambda_f=\textbf{C}:e_i\otimes e_j$, where $\textbf{C}$ is the right Cauchy-Green deformation tensor and is linked to the deformation gradient tensor $\textbf{F}$ as $\textbf{C}=\textbf{F}^T \textbf{F}$.

The fiber family was assumed to have a $\pi$ periodic von-Mises distribution: \begin{equation} \varrho(\Theta)=4\sqrt{\frac{b}{2\pi}}\frac{\displaystyle{e^{b\cos(2\Theta)+b}}}{\text{erfi}(\sqrt{2b})} \end{equation}

The characteristic invariant is defined as

\begin{align} \Lambda_2&=\frac{1}{4\pi}\int_{\Theta=0}^\pi \int_{\Phi=0}^{2\pi} \rho(\Theta) (\lambda_f^2-1)^2 d\Phi \sin(\Theta) d\Theta\\ &=\frac{1}{4\pi}\int_{\Theta=0}^\pi \varrho(\Theta) \sin(\Theta) \left[\int_{\Phi=0}^{2\pi} (\lambda_f^2-1)^2 d\Phi \right] d\Theta \end{align}

The following cell calculates the integral over $\Phi$ \begin{equation} g=\frac{1}{4\pi}\int_{\Phi=0}^{2\pi} (\lambda_f^2-1)^2 d\Phi \end{equation}





In [1]:



    


def ef(theta, phi):
    # This returns the fiber orientation vector
    a=sin(theta)
    v=vector([cos(theta), a*cos(phi), a*sin(phi)])
    return v

def lambda_f(C, theta_f, phi_f):
    # This function returns fiber stretch $\lambda_f^2$ after deformation    
        
    a=ef(theta_f, phi_f)*C*ef(theta_f, phi_f)

    return a
    

def rhom(theta, phi, b):
    # This function returns the pi periodic von Mises distribution without the normalization factor

    rho=exp(b*(cos(2*theta)+1))        
    return rho

def KKK(b):
    # This function return the normalization factor for the von Mises distribution function
    import scipy.special
    erfi=-I*scipy.special.erf(complex(I*sqrt(2*b)))

    KKK=1./(4*sqrt(b/2/pi)/erfi)

    return KKK
    

from sage.symbolic.integration.integral import definite_integral

# Define right Cauchy-Green deformation tensor
var('C11, C12, C13, C21, C22, C23, C31, C32, C33')
C=matrix(SR, 3, 3, [[C11, C12, C13], [C12, C22, C23], [C13, C23, C33]])

# $\frac{1}{4\pi} (\lambda_f^2-1)^2
var('theta_f, phi_f')
L_phi(theta_f, phi_f,b,K)=((lambda_f(C, theta_f=theta_f, phi_f=phi_f))-1)^2/4/pi

# Integral over Phi from 0 to 2\pi
g(theta_f)=definite_integral(L_phi(theta_f, phi_f), phi_f, 0, 2*pi).simplify_trig()

show(g)

1

The following cell reorganizes expression $g$ into Eqn. (18).





In [2]:



    


def reorg(func, terms):
    tmpfunc=copy(func).simplify_full()
    coefs=[]
    for each in terms:
        coef=tmpfunc.coefficient(each)
        coefs.append(coef)
        tmpfunc=(tmpfunc-coef*each).simplify_full()
    coefs.append(tmpfunc)
    return coefs


Terms=[cos(theta_f)^4, sin(theta_f)^4, cos(theta_f)^2*sin(theta_f)^2, cos(theta_f)^2, sin(theta_f)^2]
M=reorg(g, Terms)
show(M)

2

The following cell defines $\kappa_i$ as a function of $b$ and then creates Fig. 3 (Fig3.png).

The function kappas(b, terms, NI) (NI=0,4) correpsonds to $\kappa_i$ (i=1,5) in Eqn. (21).





In [3]:



    


import matplotlib.pyplot as plt
from numpy import array
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
minorLocator   = MultipleLocator(1)
from matplotlib import colors
from cycler import cycler
import matplotlib

# \kappas _ 1 to 5 (NI=0 to 4)
def kappas(b,terms,NI):
    intgrnd=rhom(theta_f, phi_f, b)*terms[NI]*sin(theta_f)    
    return numerical_integral(intgrnd, 0, pi, params=[b])[0]/KKK(b)

FF=[]
for kI in xrange(5):
    FF.append([])
    for ni in xrange(11):
        b=max(ni,1e-6)
        FF[kI].append((b, kappas(b, Terms, kI)))

        
matplotlib.rcParams.update({'font.size': 15})

plt.rc('lines', linewidth=2)

fig, ax1=plt.subplots()

from scipy import interpolate
import numpy as np

x1,y1=zip(*FF[0])
xnew = np.linspace(np.array(x1).min(),np.array(x1).max(),300)
tck = interpolate.splrep(np.array(x1), np.array(y1), s=0)
power_smooth = interpolate.splev(xnew, tck, der=0)
k1,=ax1.plot(xnew,power_smooth, 'r-', label='$\kappa_1$')

x1,y1=zip(*FF[1])
xnew = np.linspace(np.array(x1).min(),np.array(x1).max(),300)
tck = interpolate.splrep(np.array(x1), np.array(y1), s=0)
power_smooth = interpolate.splev(xnew, tck, der=0)
k2,=ax1.plot(xnew,power_smooth, 'g-.', label='$\kappa_2$', dashes=(8,2,1,2))

x1,y1=zip(*FF[2])
xnew = np.linspace(np.array(x1).min(),np.array(x1).max(),300)
tck = interpolate.splrep(np.array(x1), np.array(y1), s=0)
power_smooth = interpolate.splev(xnew, tck, der=0)
k3,=ax1.plot(xnew,power_smooth, 'b--', label='$\kappa_3$', dashes=(10,5))

x1,y1=zip(*FF[3])
xnew = np.linspace(np.array(x1).min(),np.array(x1).max(),300)
tck = interpolate.splrep(np.array(x1), np.array(y1), s=0)
power_smooth = interpolate.splev(xnew, tck, der=0)
k4,=ax1.plot(xnew,power_smooth, 'c:', label='$\kappa_4$')

x1,y1=zip(*FF[4])
xnew = np.linspace(np.array(x1).min(),np.array(x1).max(),300)
tck = interpolate.splrep(np.array(x1), np.array(y1), s=0)
power_smooth = interpolate.splev(xnew, tck, der=0)
k5,=ax1.plot(xnew,power_smooth, 'm-.', label='$\kappa_5$', dashes=(20,2,1,2))

ax1.legend([k1, k2, k3, k4, k5], ['$\kappa_1$','$\kappa_2$','$\kappa_3$','$\kappa_4$','$\kappa_5$'], bbox_to_anchor=(0.95, 0.7))
ax1.set_ylabel('$\kappa_{i\,(i=1,5)}$', fontsize=25)
ax1.set_xlabel('$b$', fontsize=25)
ax1.xaxis.set_minor_locator(MultipleLocator(1))
ax1.yaxis.set_minor_locator(MultipleLocator(2))

plt.savefig('Fig3.png')

3

Define shearing and stretching motion and calculate the corresponding right Cauchy-Green deformation tensor





In [4]:



    


var('gamma, l1')
F_Shearing=matrix(SR, 3, 3, [1, gamma, 0, 0, sqrt(1-gamma^2), 0, 0, 0, 1/sqrt(1-gamma^2)])
F_Stretch=matrix(SR, 3, 3, [l1, 0, 0, 0, 1/sqrt(l1), 0, 0, 0, 1/sqrt(l1)])
F_Total=F_Shearing*F_Stretch
C_Total=F_Total.transpose()*F_Total

4

This cell creates the comparison of strain energy function between direction integration, GST and characteristic invariant approaches, as seen in Figs. (4)-(6). The results are saved in files 'SEL???b???k2??.png'.





In [ ]:



    


def Lambda(gamma,l1,b_in, C_in):
    """
    Calculates the characteristic invariant \Lambda_2 using Eqn. (20)"""
    lmd=1
    for ni in xrange(5):
        lmd+=M[ni](C11=C_in[0][0], C12=C_in[0][1], C13=C_in[0][2],
                   C21=C_in[1][0], C22=C_in[1][1], C23=C_in[1][2],
                   C31=C_in[2][0], C32=C_in[2][1], C33=C_in[2][2])*kappas(b_in, Terms, NI=ni)
        
    return lmd

def W_CH(l1_in, b_in, k1, k2):
    r"""
    Calculates the strain energy function using the characteristic invariant"""
    W=[]
    for ni in range(21):
        gamma_in=0.01*ni
        lambda_2 = Lambda(gamma_in, l1_in, b_in, C_Total(gamma=gamma_in, l1=l1_in))
        W_ch=k1/(2*k2)*(exp(k2*lambda_2)-1)
        W.append([gamma_in,W_ch])
    return W


def strain_energy_DI(b_in, C_in, k1, k2):
    r"""
    This function calculates the strain energy function of dispersed fiber family using the direct integration approach"""
    import scipy.integrate

    wch_integral(theta_f,phi_f)=((exp(k2*(lambda_f(C=C_in, theta_f=theta_f, phi_f=phi_f)-1)^2)-1)
                                 *rhom(theta_f, phi_f,b=b_in))*sin(theta_f)/4/pi*k1/k2/2
    W=scipy.integrate.dblquad(wch_integral, 0, 2*pi, lambda x:0, lambda x:pi, epsabs=1e-5, epsrel=1e-5)

    return W[0]/KKK(b=b_in)


def W_DI(l1_in, b_in, k1, k2):
    W=[]
    for ni in range(21):
        gamma_in = 0.01*ni
        W_di=strain_energy_DI(b_in, C_Total(gamma=gamma_in, l1=l1_in), k1, k2)
        W.append([gamma_in, W_di])
    return W



def W_GST(l1_in, b_in, k1, k2):
    r"""
    Calculates the strain energy function using the GST invariant
    """
    W=[]
    import scipy.integrate
    for ni in range(21):
        gamma_in=0.01 * ni
        lambda_f_gst_itg(x,y)=((lambda_f(C=C_Total(gamma=gamma_in, l1=l1_in), theta_f=x, phi_f=y)-1))*rhom(x, y,b=b_in)*sin(x)/4/pi         
        E_i = scipy.integrate.dblquad(lambda_f_gst_itg, 0, 2*pi, lambda x:0, lambda x:pi, epsabs=1e-6, epsrel=1e-6)[0]/ KKK(b=b_in)
        W_g=k1/(2*k2)*(exp(k2*E_i^2)-1)
        W.append([gamma_in, W_g])
    return W

def Strain_Energy_Comparison(lll, bbb, k1, k2, Filename):
    ppp=[]
    labels=['CH', 'DI', 'GST']
    psi_ch=W_CH(l1_in=lll,b_in=bbb, k1=k1, k2=k2)
    ppp.append(list_plot(psi_ch, color=rainbow(5)[0], plotjoined=True, linestyle="--", marker="o", legend_label=labels[0]))

    psi_di=W_DI(l1_in=lll, b_in=bbb, k1=k1, k2=k2)
    ppp.append(list_plot(psi_di, color=rainbow(5)[1], plotjoined=True, legend_label=labels[1]))

    psi_gst=W_GST(l1_in=lll, b_in=bbb, k1=k1, k2=k2)
    ppp.append(list_plot(psi_gst, color=rainbow(5)[2], plotjoined=True, legend_label=labels[2]))

    bb=plot(sum(ppp))
    bb.save(Filename,legend_loc="upper center")
    
    

import os.path, os

LL=[1.0, 1.1, 1.2]
Bs=[1., 5., 10.]
K2s=[1, 10, 50]

for eachl in LL:
    for eachb in Bs:
        for eachk2 in K2s:
            filename='SE_L_%.3s_b_%.4s_k_%.3s.png' %(eachl, eachb, eachk2)
            Strain_Energy_Comparison(lll=eachl, bbb=eachb, k1=996, k2=eachk2, Filename=filename)

Content source: lianggecm/characteristicinvariant

Similar notebooks:

notebook.community | gallery | about