import matplotlib.pyplot as plt import numpy as np E=200e9 #Young's modulus rho=7800. #density L=1. #total length b=0.010 #beam depth h=0.005 #beam height A=b*h #beam cross sectional area I=b*h**3/12 #beam moment of inertia N=40 #number of elements Nnodes=N+1 #number of nodes ell=L/N #length of beam element Nmodes=12 #number of modes to keep mass=rho*L*A #total mass Irot=1/12*mass*L**2 #rotary inertia x=np.arange(0,ell+L,ell) K_e=E*I/ell**3*np.array([[12, 6*ell, -12, 6*ell], [6*ell, 4*ell**2, -6*ell, 2*ell**2], [-12, -6*ell, 12, -6*ell], #x2 [6*ell, 2*ell**2, -6*ell, 4*ell**2]]) #theta2 bcs = 'clamped-clamped' # 'clamped-clamped' # 'clamped-sliding' # simplysupported # 'cantilever' Input boundary conditions here loadcase = 'uniformpres' #'pointforce' # 'pointmoment' # 'unifromdistribmoment' Input load case here # 'uniformpres' # 'concmomenteachnode' K=np.zeros((2*(N+1),2*(N+1))) # M=zeros(2*(N+1),2*(N+1)) for i in range(1,N+1): K[(2*(i-1)+1-1):(2*(i-1)+4), (2*(i-1)+1-1):(2*(i-1)+4)]=K[(2*(i-1)+1-1):(2*(i-1)+4)][:,(2*(i-1)+1-1):(2*(i-1)+4)]+K_e f=np.zeros((2*Nnodes,1)) if loadcase == 'pointforce': print('Load: concentrated force at x=L/2') #load: concentrated force on center of beam f[Nnodes-1]=1 elif loadcase == 'pointmoment': print('Load: concentrated moment at x=L/2') #load: concentrated moment on center of beam f[Nnodes]=1 elif loadcase == 'uniformdistribmoment': print('Load: uniformly distributed moment') #load: uniform distributed moment on entire beam m=1.0 m_el=np.zeros((4,1)) m_el[:,0] = np.array([-m,0,m,0]) for i in range(1,N+1): f[2*(i-1)+1-1:2*(i-1)+4]=f[2*(i-1)+1-1:2*(i-1)+4]+m_el elif loadcase == 'concmomenteachnode': print('Load: Concentrated moment on each node') #load: concentrated moment on each node f[1:2*Nnodes:2]=1 elif loadcase == 'uniformpres': print('Uniform pressure') #load: uniform distributed load on entire beam q=-rho*A*9.81 q_el=np.zeros((4,1)) q_el[:,0] = np.array([q*ell/2, q*ell**2/12, q*ell/2, -q*ell**2/12]) for i in range(1,N+1): f[2*(i-1)+1-1:2*(i-1)+4]=f[2*(i-1)+1-1:2*(i-1)+4]+q_el else: print('Unknown loading') if bcs == 'clamped-clamped': print('BCs: clamped-clamped') #clamped-clamped beam BCs K=np.delete(K,[2*Nnodes-2,2*Nnodes-1],0) #K[1:2,:]=[] K=np.delete(K,[2*Nnodes-2,2*Nnodes-1],1) #K[:,1:2]=[] f=np.delete(f,[2*Nnodes-2,2*Nnodes-1]) #f[1:2]=[] K=np.delete(K,[0,1],0) #K[1:2,:]=[] K=np.delete(K,[0,1],1) #K[:,1:2]=[] f=np.delete(f,[0,1]) #f[1:2]=[] elif bcs == 'simplysupported': print('BCs: simply supported') #simply supported beam BCs K=np.delete(K,[2*Nnodes-2],0) #K[1:2,:]=[] K=np.delete(K,[2*Nnodes-2],1) #K[:,1:2]=[] f=np.delete(f,[2*Nnodes-2]) #f[1:2]=[] K=np.delete(K,[0],0) #K[1:2,:]=[] K=np.delete(K,[0],1) #K[1:2,:]=[] f=np.delete(f,[0]) #f[1:2]=[] elif bcs == 'cantilever': print('BCs: cantilever') #cantilever beam BCs K=np.delete(K,[0,1],0) #K[1:2,:]=[] K=np.delete(K,[0,1],1) #K[:,1:2]=[] f=np.delete(f,[0,1]) #f[1:2]=[] elif bcs == 'clamped-sliding': print('BCs: clamped-sliding') #clamped-sliding beam BCs K=np.delete(K,[2*Nnodes-1],0) #K[1:2,:]=[] K=np.delete(K,[2*Nnodes-1],1) #K[:,1:2]=[] f=np.delete(f,[2*Nnodes-1]) #f[1:2]=[] K=np.delete(K,[0,1],0) #K[1:2,:]=[] K=np.delete(K,[0,1],1) #K[:,1:2]=[] f=np.delete(f,[0,1]) #f[1:2]=[] else: print('Unknown boundary conditions.') dx_vec=np.linalg.solve(K,f) if bcs == 'clamped-clamped': print('BCs output: clamped-clamped') #clamped-clamped dx=np.hstack([0., dx_vec[0:2*Nnodes-5:2], 0.]) dtheta=np.hstack([0., dx_vec[1:2*Nnodes-4:2], 0.]) elif bcs == 'simplysupported': print('BCs output: simply-supported') #simply-supported dx=np.hstack([0., dx_vec[1:2*Nnodes-4:2], 0.]) dtheta=np.hstack([dx_vec[0:2*Nnodes-3:2], dx_vec[2*Nnodes-3]]) elif bcs == 'cantilever': print('BCs output: cantilever') #cantilever dx=np.hstack([0., dx_vec[0:2*Nnodes-2:2]]) dtheta=np.hstack([0., dx_vec[1:2*Nnodes-1:2]]) elif bcs == 'clamped-sliding': print('BCs output: clamped-sliding') #clamped-sliding beam BCs dx=np.hstack([0., dx_vec[0:2*Nnodes-3:2]]) dtheta=np.hstack([0., dx_vec[1:2*Nnodes-4:2], 0.]) else: print('Output: Unknown boundary conditions.') plt.figure(1) plt.subplot(211) plt.plot(x,dx) plt.ylabel('displacement [m]') #plt.title(['BCs: ' bcs 'Load case: ' loadcase]) plt.show() plt.subplot(212) plt.plot(x,dtheta) plt.ylabel('slope [radians]') plt.xlabel('X [m]')]]>

Heat transfer by radiation is nominally given by

is the heat flux (W/m in SI units), is the Stefan-Boltzmann constant (5.670 x 10 W/mK), and is the emissivity, and and are the *absolute temperatures* of the surfaces. Note: Another consideration that one must also include is the view factor which describes how well the two bodies are facing each other, but we will discuss that in a future post.

When the temperatures and are not too far apart, then we can linearize the nonlinear relation above to obtain a linear convection relation of the form

where is the convection coefficient (W/mK in SI units) to obtain

Therefore, the effective heat transfer coefficient is

where is a characteristic temperature. As the characteristic temperature, one may use the nominal value , , or some value in between, such as the mean. Depending on the problem one may choose the value based in order to make a conservative approximation.

This relation can be derived in two ways. First, one can factor the original relation into

If one defines that mean temperature as and assumes that T and T_0 are not too far apart, then and one finds

A second way involves using a Taylor series expansion with . Then the heat flux is given by

Expanding and dropping higher order terms of and above, we obtain

Above, I wrote that a few different characteristic temperatures can be used. Those can be obtained here depending on the expansion. One could write and for example.

Suppose we have two bodies with emissivities and at temperatures and transferring heat between each other as sketched below. Each body emits radiation as but *the net heat transfer is not the simple difference* because the incoming and outgoing radiation between the two bodies are coupled.* *The net heat transfer is given by

Therefore, the effective emissivity is

which can be used with the relation

One can see that if either or or both, then the more familiar relation for heat transfer is obtained.

One can find a full derivation of this relation via this link.

]]>

]]>

A great reference for its elastic properties including when isotropic approximations are appropriate is the paper by Hopcroft et al (M. A. Hopcroft, W. D. Nix, and T. W. Kenny, “What is the Young’s Modulus of Silicon?,”

Journal of Microelectromechanical Systems, vol. 19, pp. 229-238, 2010. link to researchgate https://www.researchgate.net/publication/224124507_What_is_the_Young%27s_Modulus_of_Silicon). Additionally Hopcroft’s course notes are useful http://micromachine.stanford.edu/~hopcroft/Publications/Hopcroft_E_Si_v1p1.pdf

]]>

We begin with a very brief refresher of Euler-Bernoulli beam theory as described in one of my favorite texts \emph{Mechanical Vibration} by S.S. Rao (link to book). The governing equation for a vibrating beam is

For free vibration, and separation of variables is used to obtain two equations

and

where

The mode shape of the beam is

The boundary conditions of the beam provide four equations which are used to solve for the . Fixed boundary conditions result in and . Free boundary conditions result in zero moment and shear force at the ends and . Notice here that the free boundary conditions are two additional derivatives of the fixed boundary conditions. Pinned boundary conditions have zero displacement and zero moment at the ends and . Sliding restrained boundary conditions have zero slope and zero shear force at the ends , . Again these two boundary conditions are related by two differentiations.

For the fixed-fixed beam, substitution of the boundary conditions results in a four equations

Similarly for the free-free beam

Comparing the two equations, the only differences are that the first two columns of the free-free case are the negative of the first two columns of the fixed-fixed case. To solve for the non-trivial solution, one must set the determinant of the matrix equal to zero to obtain the characteristic equation. One property of the determinant is that . When we compare the transpose of the matrices we see that the first two rows of the free-free case and the fixed-fixed case differ only by a negative sign and can be made identical. Therefore, the characteristic equation will have the same solution for in each case and the natural frequencies are identical. When one goes back to solve for the mode shapes, it is quite clear that the negative sign results in different shapes.

The result of the characteristic equation is

The mode shape for the fixed-fixed beam is

where

The mode shape for the free-free beam is

where is the same as for the fixed-fixed case given above.

Below we plot the first mode shape as well as its derivatives for the fixed-fixed and free-free boundary conditions. The maximum displacement has been normalized to 1 and and and have also been set to 1 without loss of generality.

Another useful insight is to compare the strain and kinetic energy densities. Because the natural frequencies are the same, we know that the ratio of potential energy and strain energy are the same, but the energy density shows how the beams are different.

A similar procedure can also be followed to compare pinned and sliding boundary conditions. Furthermore, one also finds that this holds true for natural frequencies of plates.

]]>

The shoelace formula gets its name from the arrangement of the coordinates and how they are combined to calculate the area. Arrange the *x-y* coordinates of the polygon in a *(n+1)x2* matrix where the order is determined by a counterclockwise pattern around the perimeter and the starting point is also repeated as the last row in the matrix.

Again, notice that the order is counterclockwise and that the first point is repeated in the last line of the matrix (there are *n+1* rows if the polygon has *n* points).

To calculate the area, we add when multiplying down and to the right and subtract pairs when multiplying down and to the left. That is

In the example polygon shown above, we have

Performing the calculation, we obtain

]]>

Several command line options for ANSYS are hard-coded into this script. To generate the appropriate options and paths for your machine, set up your solve using the ANSYS launcher and then display the command that it would invoke from “Tools –> Display Command Line.”

I don’t show it here, but one can change the arguments to the APDL script by having the Python program edit either the body of the APDL script or its input files.

from subprocess import call import datetime import os def runAPDL(ansyscall,numprocessors,workingdir,scriptFilename): """ runs the APDL script: scriptFilename.inp located in the folder: workingdir using APDL executable invoked by: ansyscall using the number of processors in: numprocessors returns the number of Ansys errors encountered in the run """ inputFile = os.path.join(workingdir, scriptFilename+".inp") # make the output file be the input file plus timestamp outputFile = os.path.join(workingdir, scriptFilename+ '{:%Y%m%d%H%M%S}'.format(datetime.datetime.now())+ ".out") # keep the standard ansys jobname jobname = "file" callString = ("\"{}\" -p ansys -dis -mpi INTELMPI" " -np {} -dir \"{}\" -j \"{}\" -s read" " -b -i \"{}\" -o \"{}\"").format( ansyscall, numprocessors, workingdir, jobname, inputFile, outputFile) print("invoking ansys with: ",callString) call(callString,shell=False) # check output file for errors print("checking for errors") numerrors = "undetermined" try: searchfile = open(outputFile, "r") except: print("could not open",outputFile) else: for line in searchfile: if "NUMBER OF ERROR" in line: print(line) numerrors = int(line.split()[-1]) searchfile.close() return(numerrors) def main(): ansyscall = "C:\\Program Files\\ANSYS Inc\\v180\\ansys\\bin\\winx64\\MAPDL.exe" numprocessors = 8 workingdir = "G:\\scriptAnsDir" scriptFilename = "dymmyAnsysScript" nErr = runAPDL(ansyscall, numprocessors, workingdir, scriptFilename) print ("number of Ansys errors: ",nErr) if __name__ == '__main__': main()]]>

One topic that is particularly interesting and well written is the section on subcritical crack growth. Subcritical crack growth (also called stress corrosion) is the growth or extension of a crack over time with stress intensity factor* less than* the critical stress intensity factor . This is interesting because in classic linear elastic fracture mechanics (LEFM) cracks in brittle materials are typically viewed as stable (no growth under constant load) if the the stress intensity factor is less than the the critical value . The text only briefly mentions the underlying cause of subcritical crack growth in terms of chemical bonds breaking in the neighborhood of the crack tip.

The crack growth is given by

where is the crack growth velocity with respect to time, is the length of the crack, is time, is the stress intensity factor, is the critical stress intensity factor, and and are constants based on the particular material, temperature and environmental conditions. Another common parameter is which is given by

where is a constant related to the geometry, typically of order 1, with value for and infinite plate.

Zerodur (manufacturer Schott’s website) is one of the most commonly used ceramics or glass ceramics by precision engineers. Subcritical crack growth is an important consideration for Zerodur in the presence of water or water vapor. The term stress corrosion is often used in this case as many like to think of the water molecules attacking the bonds at the root of the crack.

“Fracture Toughness and Crack Growth of Zerodur” (*NASA Technical Memorandum 4185*) by Michael J. Viens, April 1990 (link to document) gives the stress corrosion constants to be and with 100% humidity. *The Properties of Optical Glass* edited by Hans Bach and Norbert Neuroth (google books link, amazon.com link) gives the value of for Zerodur. They also give the values of for Duran and soda-lime glass as 22.4 and 18.1, respectively. The Zerodur technical document by Schott TIE-33 (link) has an extensive discussion, examples, and list of fracture and related properties including stress corrosion.

]]>

]]>

Consider as shown in the figure above a c-shaped iron core with a coil wrapped around one side. The iron core has a small air gap, into which a copper plate is inserted. If we drive a constant current through the coil, we induce a magnetic flux in the iron. The DC flux density across the air gap in our example is about 0.15 T. For a static (DC) field, the copper behaves just like air, not altering the magnetic field.

But when the current and therefore the magnetic field vary harmonically with time, an electromotive force is induced according to Faraday’s Law

(1)

where is the induced electric field, is an element of the boundary of the area through which the magnetic flux passes. We can obtain by integrating the flux density over the area according to

(2)

For our example, the flux density is nearly uniform over the area and drops off rapidly outside of , and hence is well approximated by . If varies harmonically with frequency , then the magnitude of and hence the magnitude of the EMF around any loop surrounding the area is .

The current density induced by the EMF on a circular loop of radius in a conductor is given by , where is the resistivity of the conductor. So the magnitude of the current density is approximately

(3)

Thus, for our example with =0.15 T and =(5 mm) and =1.7e-8 -m, we expect the flux density just outside the 5 mm square (at =3 mm) to be 2.8e6 A/m at 40 Hz. This agrees pretty well with the numerical results shown in the animation above, which was generated using Ansoft Maxwell from ANSYS.

This simple approximation works quite well at low frequencies. At higher frequencies, the EMF of the induced currents prevents the magnetic field from penetrating the conductor. The currents drop off rapidly with a skin depth characterized by

(4)

where is the magnetic permeability of the material. For copper at 40~Hz, we have a skin depth of about 10 mm, much larger then the 1 mm thickness of the plate in our example. We can therefore safely consider the induced current to be uniform through the thickness of the plate.

]]>