Wolfram Language for Materials Sciences applications

Gaurav Vishal, Consultant, WTG

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

OUTLINES

Curated Physical, Chemical and other data of Elements and Molecules

ElementData

ChemicalData

Molecule

LatticeData

Automatic Detection of Image Features

Calculating Size distribution of nanoparticles

Calculating Phase fraction in Microstructure

Solving PDEs

Diffusion Problem

Heat Transfer Problem

Finite Element Programming (Modelling Phase Change in Materials )

Element Data in Wolfram Language

Element Data in Wolfram Language

ElementData provides phyisical, chemical and other properties of 118 elements which are present in periodic table

ElementData[]

ElementData[]//Length

ElementData["Properties"]

ElementData["O","AtomicWeight"]

Precision[%]

ElementData["Tungsten","ThermalConductivity"]

ElementData[Entity["Element",#],"Electronegativity"]&/@{"F","Cl","Br","I"}

ElementData["Hydrogen","StableIsotopes"]

ElementData["Copper","ElectronConfiguration"]

ElementData["Copper","CrystalStructure"]

%["Image"]

Chemical Data in Wolfram Language

Chemical Data in Wolfram Language

ChemicalData provides features curated data for over 80k molecules, drawn from authoritative sources such as NIST, PubChem, ChemSpider

ChemicalData["Properties"]

ChemicalData["Properties"]//Length

ChemicalData["Benzene","CombustionHeat"]

ChemicalData["Benzene","AtomPositions"]

ChemicalData["Benzene","MolarMass"]

ChemicalData["Benzene","AlternateStandardNames"]

ChemicalData[#,"IUPACName"]&/@%

Molecule in Wolfram Language

Molecule in Wolfram Language

Gives the molecule entity corresponding to the input name

Borazine (Inorganic Benzene)

mol1=Molecule["borazine"]

Table[MoleculePlot[mol1,PlotTheme->theme,],{theme,{"AllAtom","HeavyAtom","Aromatic","Monochrome"}}]

AtomList[mol1]

BondList[mol1]

Create a molecule from a systematic chemical name

mol2=Molecule["N,N-dimethyl-3-[(9R,10S)-10-methyl-2-(trifluoromethyl)-9,10-dihydroanthracen-9-yl]propan-1-amine"]

MoleculePlot[mol2]

Highlight all carbon-Florine bonds:

MoleculePlot[mol2,Bond[{"C","F"}]]

Create a molecule from Atom and Bond list

mol3=Molecule[{"O","H","H"},{Bond[{1,2},"Single"],Bond[{1,3},"Single"]}]

Replace Atom O with S

MoleculeModify[mol3,{"ReplaceAtom",1->"S"}]

Molecule properties

MoleculeValue["Properties"]

mol3["MolecularMass"]

MoleculeValue[mol3,"GeometricProperties"]

Stereochemistry Elements

Create a molecule with a defined stereo-center:

m1 = Molecule[{Atom["N"],Atom["C"],Atom["Br"],Atom["Cl"],Atom["C"]}, {Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{2,4},"Single"],Bond[{2,5},"Single"]},StereochemistryElements->{\[LeftAssociation]"StereoType"->"Tetrahedral","ChiralCenter"->2,"Direction"->"Counterclockwise","FiducialAtom"->1,"Ligands"->{3,4,5}\[RightAssociation]}]

MoleculePlot3D[m1]

m1[{"AtomChirality",2}]

Create its enantiomer:

m2 = Molecule[{"N","C","Br","Cl","C"},{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{2,4},"Single"],Bond[{2,5},"Single"]},StereochemistryElements->{\[LeftAssociation]"StereoType"->"Tetrahedral","ChiralCenter"->2,"Direction"->"Clockwise","FiducialAtom"->1,"Ligands"->{3,4,5}\[RightAssociation]}]

MoleculePlot3D[m2]

m2[{"AtomChirality",2}]

Molecule Recognize

img=;

ImageQ[img]

recognizedMol=MoleculeRecognize[img]

MoleculeValue[recognizedMol,"MolarMass"]

MoleculePlot[recognizedMol,PlotTheme->"Aromatic"]

Import/ Export formats

SDF - Structure-data File

XYZ import - connecting the dots

CIF - Crystallographic Information File

HIN - HyperChem Format

CML - Chemical Markup Language

FCHK - Formatted Checkpoint File

CUBE - Volumetric Orbital Data

cml=Import["D:\\Projects\\Material Science\\Show & Tell\\toluene.cml",{"Molecule",1}]

MoleculePlot[cml]

cml=Import["D:\\Projects\\Material Science\\Show & Tell\\toluene.cml",{"Metadata",1}]

LatticeData in Wolfram Language

LatticeData

LatticeData[]

LatticeData["Properties"]

LatticeData["FaceCenteredCubic","Image"]

LatticeData["FaceCenteredCubic","AlternateNames"]

LatticeData["FaceCenteredCubic","PackingRadius"]

LatticeData["FaceCenteredCubic","Density"]

Calculating Size distribution of nanoparticles

Size distribution of nanoparticles

A transmission electron microscopy (TEM) image shows the iron oxide nanoparticles coated in a mesoporous silica that are used in the tissue warming process.

temImg=;

Binarized Image

binarizedImg=Binarize[temImg]

Color Negate

colorNeg=ColorNegate[binarizedImg]

MorphologicalComponents

pixelArray=MorphologicalComponents[colorNeg];

Delete Particles on border

cleanedPixelArray=DeleteBorderComponents[pixelArray];

Generate Color Image from cleanedPixelArray

Colorize[cleanedPixelArray]

Delete small particles

cleanedPixelArray1=DeleteSmallComponents[cleanedPixelArray,25];
Colorize[cleanedPixelArray1]

Measurements

equivDiskRadii=ComponentMeasurements[cleanedPixelArray1,"EquivalentDiskRadius"][[All,2]]

Note that the values are in pixels, so we would need to calculate the

equivDiskRadiiNanometers=equivDiskRadii*20/60

Calculating  pixel size using scale bar

temImg

whitePixel=PixelValuePositions[temImg,White];
selectedWhitePixel=Select[whitePixel,SameQ[Last[#], Min[whitePixel[[All,2]]]]&];

HighlightImage[temImg,selectedWhitePixel,Method->{"DiskMarkers",10}]

Length[selectedWhitePixel]

100  nm scale bar is about ≈ 40 pixels in length, converting the equivalent disc radii to nanometers gives:

equivDiskRadiiNanometers=equivDiskRadii*100/40

Particle Size Distribution:

Histogram[equivDiskRadiiNanometers]

Calculating Phase fraction in Microstructure

Calculating Phase fraction in Microstructure

Microstructure of Hypo-eutectic white cast (iron^1)

whiteCastIron=;

Binarized Image

binarizedMicrograph=Binarize[whiteCastIron]

ColorNegate

colorNegatedMicrograph=ColorNegate[binarizedMicrograph]

Measurements

Using Pixel Count

pixelValues=Flatten[ImageData[colorNegatedMicrograph]]

pixelTally=Tally[pixelValues]

totalPixel=Total[pixelTally[[All,2]]]

percentCementite=N[((pixelTally[[2,2]])/(Total[pixelTally[[All,2]]]))]

Using ComponentMeasurements

pixelArray=MorphologicalComponents[colorNegatedMicrograph];
cleanedPixelArray=DeleteSmallComponents[pixelArray,25];

Colorize[cleanedPixelArray]

equivArea=ComponentMeasurements[cleanedPixelArray,"Area"];
totalCementite=Total[Values[equivArea]]

percentCementite1=totalCementite/totalPixel

Using Lever rule:

% Cementite = ((6.67 - wt % C )/(6.67-0.76))

Solve[0.515716`==((6.67-wtC)/(6.67-0.76)),wtC]

Solving PDEs

PDEs

Some Categories of PDEs of interest to material scientist are parabolic equation with the form:

which is a linear second order homogeneous parabolic PDE with constant coefficient. 

These parabolic equation have a derivative that is second order with respect to variable x.

Boundary Condition

Solving parabolic PDEs require an initial condition at t = 0 and two boundary conditions.

Dirichlet boundary condition

Conditions that specify the solution at the boundary.

NDSolveValue[{(Subscript[∇, {x,y}]^2)u[x,y]==0,DirichletCondition[u[x,y]==Sin[x y],True]},u,{x,y}∈Disk[]]

Neumann boundary condition

Conditions that specify the derivative of the solution at the boundary.

NDSolveValue[{-(Subscript[∇, {x,y}]^2)u[x,y]==NeumannValue[1.,x≥0.35],DirichletCondition[u[x,y]==0.,x≤-0.3]},u,{x,y}∈Disk[]]

Solving (2^nd) order PDEs (Heat & Diffusion Equations)

Heat Equation

heqn=D[u[x,t],t]== D[u[x,t],{x,2}]

Diffusion Equation

1 Dimension Diffusion eqn

deqn=D[c[x,t],t]== D[c[x,t],{x,2}]

Diffusion with First order Chemical Reaction

deqn2=D[c[x,t],t]== D[c[x,t],{x,2}] -k*c[x,t]

Diffusion convection equation

deqn3=D[c[x,t],t]== D[c[x,t],{x,2}] -ν*D[c[x,t],{x,1}]

Spherical Diffusion

deqn4=D[c[r,t],t]== D[c[r,t],{r,2}] -(2/r)*D[c[r,t],{r,1}]

1D Heat Transfer Problem

heqn=D[u[x,t],t]==D[u[x,t],{x,2}]

Initial Condition

ic=u[x,0]==E^(-x^2);

sol=DSolveValue[{heqn,ic},u[x,t],{x,t}]

Plot[Evaluate[Table[sol,{t,0,4}]],{x,-5,5},PlotRange->All,Filling->Axis]

Initial value problem for the heat equation with piecewise initial data.

ic=u[x,0]==UnitBox[x];

sol=DSolveValue[{heqn,ic},u[x,t],{x,t}]

Plot3D[sol,{x,-2,2},{t,0,1},PlotRange->All,PlotPoints->250,Mesh->None]

3D Heat Transfer Problem

Heat transfer equation in 3D

where u is Temperature and α is thermal diffusivity

eqn=D[u[x,y,z,t],t]== Derivative[0,0,2,0][u][x,y,z,t]+Derivative[0,2,0,0][u][x,y,z,t]+Derivative[2,0,0,0][u][x,y,z,t]

Initial Condition

init=u[x,y,z,0]==304;

bc1=DirichletCondition[u[x,y,z,t]==293,x==0];
bc2=DirichletCondition[u[x,y,z,t]==293,x==50];
bc3=DirichletCondition[u[x,y,z,t]==293,y==0];
bc4=DirichletCondition[u[x,y,z,t]==293,y==180];
bc5=DirichletCondition[u[x,y,z,t]==280,z==0];
bc6=DirichletCondition[u[x,y,z,t]==293,z==40];

usol=NDSolveValue[{eqn,init,bc1,bc2,bc3,bc4,bc5,bc6},u,{x,0,50},{y,0,180},{z,0,40},{t,0,600}]

ContourPlot3D[usol[x,y,z,40],{x,0,50},{y,0,180},{z,0,40},Contours->{294}]

ContourPlot3D[usol[x,y,z,40],{x,0,50},{y,0,180},{z,0,40},RegionFunction->Function[{x,y,z},x<25||y>90||z<20],Contours->{296,298,300,302}]

For more details on solving PDE, refer to:

https://reference.wolfram.com/language/tutorial/DSolvePartialDifferentialEquations.html

1 Dimension Diffusion Problem

deqn=D[c[x,t],t]== D[c[x,t],{x,2}]

Initial Value

iv=c[x,0]==0;

Boundary Condition

bc ={c[0,t] ==Sin[t], c[1,t]==0};

s11=NDSolve[{deqn,iv,bc},c[x,t],{x,0,1},{t,0,10}]

Plot3D[c[x,t]/.s11,{x,0,1},{t,0,10},AxesLabel->Automatic]

Using periodic boundary conditions

bc={c[0,t]==Sin[t],c[2Pi,t]==Sin[t]}

s12=NDSolve[{deqn,iv,bc},c[x,t],{x,0,2Pi},{t,0,10}]

Plot3D[c[x,t]/.s12,{x,0,1},{t,0,10}]

Finite Element Programming (Modelling Phase Change)

Finite Element Programming

Low level interface that  gives extensive control for each part of the solution process

To use the finite element functions, the NDSolve`FEM` package needs to be loaded:

Needs["NDSolve`FEM`"]

The solution can be found in three stages :

Initialization

Discretization

Solving

Phase Change in Material

Phase change follows nonlinear heat equation, consider



with temperature, density, temperature - dependent heat capacity and temperature - dependent conductivity .

Transient PDEs with Nonlinear Transient Coefficients

In order to keep things simple here, lets model the phase change in a 1 D domain

Lets assume that phase change from one phase to the next make use of smooth transition function :

SmoothedStepFunction[fmin_,fmax_,ts_,m_]:=Function[t,(fmin*Exp[ts*m]+fmax*Exp[m*t])/(Exp[ts*m]+Exp[m*t])]
θ=SmoothedStepFunction[1,0,0,50];

Show[Plot[θ[T],{T,-1,1},],Graphics[]]

Set up the nonlinear coefficients and parameters:

k[T_]:=θ[T]*2+(1-θ[T])*1/2;
Cp[T_]:=θ[T]+(1-θ[T])*2;
ρ=2*10^6;

Specify a 1 D domain:

Ω=Line[{{0},{0.01}}];

BCs:

On the left boundary there is a constant inflow of energy and on the right boundary a fixed temperature is set up. Over time, the initial phase of the material will change from the left to the right end of the domain.

bcLeft=NeumannValue[5000,x==0];
bcRight=DirichletCondition[T[t,x]==-10,x==0.01];

Set up the numerical region, the variable and solution data :

Needs["NDSolve`FEM`"]

nr=ToNumericalRegion[Ω];
vd=NDSolve`VariableData[{"DependentVariables"->{T},"Space"->{x},"Time"->t}];
sd=NDSolve`SolutionData[{"Space"->nr,"Time"->0.}];

Initialize the PDE coefficient, boundary conditions and the method data

initCoeffs=InitializePDECoefficients[vd,sd,"DiffusionCoefficients"->{{-k[T[t,x]]*IdentityMatrix[1]}},"DampingCoefficients"->{{ρ*Cp[T[t,x]]}}];
initBCs=InitializeBoundaryConditions[vd,sd,{{bcLeft,bcRight}}];
methodData=InitializePDEMethodData[vd,sd,Method->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->0.01/50}}];

Compute the stationary components of the PDE coefficients and boundary conditions .

sdpde=DiscretizePDE[initCoeffs,methodData,sd,"Stationary"];
sbcs=DiscretizeBoundaryConditions[initBCs,methodData,sd,"Stationary"];

Evaluate the initial conditions such that the values are set to everywhere

init=EvaluateOnElementMesh[{x},-10,methodData["ElementMesh"]]["ValuesOnGrid"];

Create a function to compute the discretized PDE in each time step :

discretizePDEResidual[t_?NumericQ,u_?VectorQ,dudt_?VectorQ]:=
Module[{l,s,d,tdpde,tbcs,nldpde,nlbcs,sdTemp},
	NDSolve`SetSolutionDataComponent[sd,"Time",t];
	NDSolve`SetSolutionDataComponent[sd,"DependentVariables",u];
	l=sdpde["LoadVector"];
	s=sdpde["StiffnessMatrix"];
	d=sdpde["DampingMatrix"];
	tdpde=DiscretizePDE[initCoeffs,methodData,sd,"Transient"];
	tbcs=DiscretizeBoundaryConditions[initBCs,methodData,sd,"Transient"];
	l+=tdpde["LoadVector"];
	s+=tdpde["StiffnessMatrix"];
	d+=tdpde["DampingMatrix"];
	nldpde=DiscretizePDE[initCoeffs,methodData,sd,"Nonlinear"];
	nlbcs=DiscretizeBoundaryConditions[initBCs,methodData,sd,"Nonlinear"];
	l+=nldpde["LoadVector"];
	s+=nldpde["StiffnessMatrix"];
	d+=nldpde["DampingMatrix"];
	DeployBoundaryConditions[{l,s,d},nlbcs];
	DeployBoundaryConditions[{l,s,d},tbcs];
	DeployBoundaryConditions[{l,s,d},sbcs];
	d.dudt+s.u-l
]

Solve the equations

tEnd=400;
Monitor[tTfun=NDSolveValue[{discretizePDEResidual[t,T[t],T'[t]]==0,T[0]==init},T,{t,0,tEnd},Method->{"EquationSimplification"->"Residual"},EvaluationMonitor:>(monitor=Row[{"t = ",CForm[t]}])],monitor]

Create an InterpolatingFunction from the time steps

Tfun=ElementMeshInterpolation[{tTfun["Coordinates"][[1]],methodData["ElementMesh"]},Partition[tTfun["ValuesOnGrid"],1]]

Visualize the solution

nframes=30;
frames=Table[Labeled[Plot[Tfun[t,x],{x,0,0.01},PlotRange->{-20,50},PlotTheme->"Detailed",PlotLabel->Style["T(t,x)",14],FrameLabel->{"x","Temperature"},LabelStyle->Directive[30],PlotLegends->None,ImageSize->1000],Text[Style[Row[{"Time"," = ",N[t,4]}],30]]],{t,0,tEnd,tEnd/nframes}];
ListAnimate[frames,SaveDefinitions->True,AnimationRunning->False]

For more details on FiniteElementProgramming, refer to:

https://reference.wolfram.com/language/FEMDocumentation/tutorial/FiniteElementProgramming.html

References

https://www.doitpoms.ac.uk/miclib/micrograph.php?id=66

https://nanocomposix.com/pages/magnetite-nanoparticles#target