Added splitting the forces into their contribution. There seems to be a sign... (c9250ea3) · Commits · Publications / Fused Filament Modeling - Generalization and Nonisthermal Solutions

force_plots.py

+69 −37

Original line number	Diff line number	Diff line
		@@ -9,6 +9,8 @@ import matplotlib.pyplot as plt

		from forces import F_newt_isoth,F_newt_vft,F_newt_eyr

		R = 8.314

		def geometry_1():
		alpha= 30 # degree
		Ro= 0.002/2 # m
		@@ -17,19 +19,17 @@ def geometry_1():

		return alpha,Ro,Ri,Lc

		def PP_values():

		km = # m s /K
		Th = # K
		Tm = # K
		rho_s = # kg m**(-3)
		rho_m = # kg m**(-3)
		h = # J/mol (0 für PS)
		Cs = # J/ mol K
		Q = # J/mol
		eta0 = # Pas
		def PP_material():
		km = 0# m s /K
		Tm = 0# K
		rho_s = 0# kg m**(-3)
		rho_m = 0# kg m**(-3)
		h = 0# J/mol (0 für PS)
		Cs = 0# J/ mol K
		Q = 0# J/mol
		eta0 = 0# Pas

		return km,Th,Tm,rho_s,rho_m,h,Cs,Q,eta0
		return km,Tm,rho_s,rho_m,h,Cs,Q,eta0

		def PS_material():

		@@ -39,11 +39,11 @@ def PS_material():
		rho_m = 917 # kg m**(-3)
		h = 0 # J/mol
		Cs = 1630 # J/ mol K
		eps = # K
		Tv = 308.1 # K
		eta0 = 2.27 * 10**3 # Pas
		eps = 2303.99 # K
		Tv = 287.03 # K
		eta0 = 0.08549 # Pas

		return km,Th,Tm,rho_s,rho_m,h,Cs,T0,eps,Tv,eta0
		return km,Tm,rho_s,rho_m,h,Cs,eps,Tv,eta0

		def PS_settings():
		Th = 513.15 # K
		@@ -52,7 +52,7 @@ def PS_settings():
		return Th,T0

		def PP_settings():
		Th = # K
		Th = 0 # K
		T0 = 293.15 # K

		return Th,T0
		@@ -62,7 +62,9 @@ def vft_forces(Us,
		geometry,
		settings,
		debug_iso=False,
		debug_non=False):
		split_iso=False,
		debug_non=False,
		split_non=False):

		km,Tm,rho_s,rho_m,h,Cs,eps,Tv,eta0 = material()
		Th,T0 = settings()
		@@ -74,12 +76,14 @@ def vft_forces(Us,
		# calculate isothermal
		F_iso = F_newt_isoth(km=km, Th=Th, Tm=Tm, rho_s=rho_s, rho_m=rho_m, h=h,
		Cs=Cs, T0=T0, Us=Us,
		eta0=eta0*np.exp(eps/((Th+Tm)/2 -Tv)),
		eta_eff=eta0*np.exp(eps/((Th+Tm)/2 -Tv)),
		alpha=alpha, Ro=Ro, Ri=Ri, etac=etac, Lc=Lc,
		debug=debug_iso)
		debug=debug_iso,
		split=split_iso)
		F_non = F_newt_vft(km, Th, Tm, rho_s, rho_m, h, Cs, T0,
		Us, eta0, eps, Tv, alpha, Ro, Ri, etac, Lc,
		debug=debug_non)
		debug=debug_non,
		split=split_non)

		return F_iso,F_non

		@@ -91,7 +95,7 @@ def eyr_forces(Us,
		debug_non=False):

		# extract material data, geometry data and process settings
		km,Th,Tm,rho_s,rho_m,h,Cs,Q,eta0 = material()
		km,Tm,rho_s,rho_m,h,Cs,Q,eta0 = material()
		Th,T0 = settings()
		alpha,Ro,Ri,Lc = geometry()

		@@ -101,7 +105,7 @@ def eyr_forces(Us,
		# calculate isothermal
		F_iso = F_newt_isoth(km=km, Th=Th, Tm=Tm, rho_s=rho_s, rho_m=rho_m, h=h,
		Cs=Cs, T0=T0, Us=Us,
		eta0=eta0np.exp(Q/(R(Th+Tm)/2)),
		eta_eff=eta0np.exp(Q/(R(Th+Tm)/2)),
		alpha=alpha, Ro=Ro, Ri=Ri, etac=etac, Lc=Lc,
		debug=debug_iso)

		@@ -115,9 +119,6 @@ def eyr_forces(Us,

		if __name__=="__main__":


		etac = eta0 * np.exp(Q/(R*Th))

		#
		fig,axs = plt.subplots(1,2)

		@@ -128,28 +129,59 @@ if __name__=="__main__":
		F_iso,F_non = vft_forces(Us,
		PS_material,
		geometry_1,
		PS_settings)
		PS_settings,
		debug_iso=False,
		split_iso=False,
		debug_non=False,
		split_non=False)

		# plot isoth. plot
		axs[0].plot(Us10*3,F_iso,label="isoth.")
		axs[0].plot(Us10*3,F_iso,color="k",label="isoth.")

		# plot nonisoth. plot
		axs[0].plot(Us10*3,F_iso,label="nonisoth.")
		axs[0].plot(Us10*3,F_non,color="b",label="nonisoth.")

		# calculate force contributions
		F_iso,F_non = vft_forces(Us,
		PS_material,
		geometry_1,
		PS_settings,
		debug_iso=False,
		split_iso=True,
		debug_non=False,
		split_non=True)

		# plot individual force contributions
		axs[0].plot(Us10*3,F_iso[0],color="k",
		linestyle="--",label=r"$F_{p}$")
		axs[0].plot(Us10*3,F_iso[1],color="k",
		linestyle="-.",label=r"$F_{pi}$")
		axs[0].plot(Us10*3,F_iso[2],color="k",
		linestyle=":",label=r"$F_{\tau}$")
		axs[0].plot(Us10*3,F_non[0],color="b",
		linestyle="--",label=r"$F_{p}$")
		axs[0].plot(Us10*3,F_non[1],color="b",
		linestyle="-.",label=r"$F_{pi}$")
		axs[0].plot(Us10*3,F_non[2],color="b",
		linestyle=":",label=r"$F_{\tau}$")


		# calculate forces
		F_iso,F_non = eyr_forces(Us,
		PP_material,
		geometry_1,
		PP_settings)
		#F_iso,F_non = eyr_forces(Us,
		# PP_material,
		# geometry_1,
		# PP_settings)

		# plot isoth. plot
		axs[1].plot(Us10*3,F_iso,label="isoth.")
		#axs[1].plot(Us10*3,F_iso,label="isoth.")
		# plot nonisoth. plot
		axs[1].plot(Us10*3,F_iso,label="nonisoth.")
		#axs[1].plot(Us10*3,F_iso,label="nonisoth.")

		for i in range(2):
		axs[i].legend()
		axs[i].set_xlabel(r"$U_{s} (mm/s)$")
		axs[i].set_ylabel(r"$U_{F} (N)$")
		axs[i].set_yscale("log")
		axs[i].set_yscale("symlog")
		axs[0].set_title("PS")
		axs[1].set_title("PP")
		plt.show()
		No newline at end of file

forces.py

+63 −25

Original line number	Diff line number	Diff line
		@@ -11,10 +11,10 @@ from integrals import eta_1,int_eta_dz,int_etaz_dz,intint_eta_dzdz,intint_etaz_d
		def get_delta(km, Th, Tm, rho_s, Us, h, Cs, T0):
		return (km * (Th - Tm)) / (rho_s * Us * (h + Cs * (Tm - T0)))

		def Ft_newt_isoth(U0, eta0, delta, alpha, Ro, Ri):
		def Ft_newt_isoth(U0, eta_eff, delta, alpha, Ro, Ri):


		prefac = 2 * np.pi * (eta0 * U0) / delta
		prefac = 2 * np.pi * (eta_eff * U0) / delta
		term1 = 3 * ((1/np.cos(alpha)) * (delta*(-1)) \
		((2/3) * Ro*3 - Ri Ro*2 + (1/3) Ri**3)\
		+ np.tan(alpha) * (Ro*2 /2 - 2RoRi + Ro*2 /2))
		@@ -23,24 +23,27 @@ def Ft_newt_isoth(U0, eta0, delta, alpha, Ro, Ri):

		return prefac*(term1 + term2)

		def Fp_newt_isoth(U0, eta0, delta, alpha, Ro, Ri, etac, Lc,
		debug):
		def Fp_newt_isoth(U0, eta_eff, delta, alpha, Ro, Ri, etac, Lc,
		debug, split):

		#
		inlet_p = 8np.pietacU0Lc(Ro / Ri)*4

		Fp=3np.pi(U0eta0)/(delta2 np.cos(alpha)*3) \
		Fp=3np.pi(U0eta_eff)/(delta2 np.cos(alpha)*3) \
		(((2np.log(Ro/Ri) - 3/2) Ro*4 + 2Ri*2 Ro2 + Ri4 /2)/delta \
		+ np.sin(alpha) * ((2np.log(Ro/Ri) - 7/3)Ro*3 + 2RiRo*2 +\
		Ri*2Ro - 2/3 Ri*3)) + inlet_p
		Ri*2Ro - 2/3 Ri*3))
		if debug:
		print("Pi",inlet_p)

		return Fp
		if not split:
		return Fp + inlet_p
		else:
		return Fp, inlet_p

		def F_newt_isoth(km, Th, Tm, rho_s, rho_m, h, Cs, T0,
		Us, eta0, alpha, Ro, Ri, etac, Lc,
		debug=False):
		Us, eta_eff, alpha, Ro, Ri, etac, Lc,
		debug=False,
		split=False):
		"""


		@@ -64,7 +67,7 @@ def F_newt_isoth(km, Th, Tm, rho_s, rho_m, h, Cs, T0,
		room temperature.
		Us : float
		feeder velocity.
		eta0 : float
		eta_eff : float
		Newtonian viscosity.
		alpha : float
		angle.
		@@ -76,11 +79,24 @@ def F_newt_isoth(km, Th, Tm, rho_s, rho_m, h, Cs, T0,
		viscosity in capillary.
		Lc : float
		capillary length.
		debug : boolean
		if True, print out several intermediate results for debugging.
		split : boolean
		if True, return different force contributions.

		Returns
		-------
		F: float
		feeder force.
		F: np.array float
		feeder force. Only provided if split not True
		Fp: np.array float
		pressure part of the force not stemming from the end capillary.
		Only provided if split is True
		Fpi: np.array float
		pressure part of the force stemming from the end capillary.
		Only provided if split is True
		Ftau: np.array float
		deviatoric stress part of the force.
		Only provided if split is True

		"""

		@@ -99,13 +115,16 @@ def F_newt_isoth(km, Th, Tm, rho_s, rho_m, h, Cs, T0,
		"U0",U0,
		"delta",delta)

		Ft = Ft_newt_isoth(U0, eta0, delta, alpha, Ro, Ri)
		Fp = Fp_newt_isoth(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, debug)
		Ft = Ft_newt_isoth(U0, eta_eff, delta, alpha, Ro, Ri)
		Fp = Fp_newt_isoth(U0, eta_eff, delta, alpha, Ro, Ri, etac, Lc,
		debug, split)
		if debug:
		print("Ft",Ft)
		print("Fp",Fp)

		if not split:
		return np.cos(alpha) * Fp - np.sin(alpha) * Ft
		else:
		return np.cos(alpha) * Fp[0], np.cos(alpha) * Fp[1], - np.sin(alpha) * Ft

		def eta(z,eta0,A,B):
		return eta0 * np.exp(A/(B+z))
		@@ -140,7 +159,7 @@ def Ft_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, A, B):
		return Ft

		def Fp_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, A ,B,
		debug):
		debug, split):

		# inlet pressure
		inlet_p = 8np.pietacU0Lc(Ro * 4)/(Ri**4)
		@@ -148,16 +167,19 @@ def Fp_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, A ,B,
		Fp =(-np.pi/4)eta0H(delta,A,B,delta)*(-1) np.cos(alpha)**(-1) \
		(U0(Ro*4 (2np.log(Ro/Ri) - 3/2) + 2Ri*2 Ro2 - Ri4 /2) \
		- 2W(delta,U0,alpha,A,B,delta)((2np.log(Ro/Ri) - 7/3)Ro**3 + \
		2RiRo2 + Ri2 Ro - 2/3 Ri**3)) + inlet_p
		2RiRo2 + Ri2 Ro - 2/3 Ri**3))

		if debug:
		print("Pi",inlet_p)

		return Fp
		if not split:
		return Fp + inlet_p
		else:
		return Fp, inlet_p

		def F_newt_vft(km, Th, Tm, rho_s, rho_m, h, Cs, T0,
		Us, eta0, eps, Tv, alpha, Ro, Ri, etac, Lc,
		debug=False):
		debug=False, split=False):
		"""
		Force for the Vogel-Fulcher-Tamann fluid. eta = eta0 exp(eps/(T-Tv))

		@@ -197,11 +219,24 @@ def F_newt_vft(km, Th, Tm, rho_s, rho_m, h, Cs, T0,
		viscosity in capillary.
		Lc : float
		capillary length.
		debug : boolean
		if True, print out several intermediate results for debugging.
		split : boolean
		if True, return different force contributions.

		Returns
		-------
		F: float
		feeder force.
		Fp: np.array float
		pressure part of the force not stemming from the end capillary.
		Only provided if split is True
		Fpi: np.array float
		pressure part of the force stemming from the end capillary.
		Only provided if split is True
		Ftau: np.array float
		deviatoric stress part of the force.
		Only provided if split is True

		"""

		@@ -226,13 +261,16 @@ def F_newt_vft(km, Th, Tm, rho_s, rho_m, h, Cs, T0,
		"B",B)
		Ft = Ft_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, A, B)

		Fp = Fp_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, A, B, debug)
		Fp = Fp_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, A, B,
		debug, split)

		if debug:
		print("Ft",Ft,"\n",
		"Fp",Fp)

		if not split:
		return np.cos(alpha) * Fp - np.sin(alpha) * Ft
		else:
		return np.cos(alpha) * Fp[0], np.cos(alpha) * Fp[1], - np.sin(alpha) * Ft

		def F_newt_eyr(km, Th, Tm, rho_s, rho_m, h, Cs, T0,
		Us, eta0, Q, R, alpha, Ro, Ri, etac, Lc,