This code uses an estimate of dissonance based on frequency difference between pair of partials, and determines a net "dissonance penalty" by summing over all pairs that are within some range of each other, applying the penalty for each pair.  The penalty function is key here.  At first we just used the absolute distance in frequency,  as an argument to the penalty function.    The fall off of the dissonance penalty with frequency difference of the partials should be more gradual at higher frequencies, approximately proportional to frequency above 1000 Hz, (see http://en.wikipedia.org/wiki/File:ERB_vs_frequency.svg) The  critical bandwidth is larger, expressed as frequency distance.

(*Load these packages*)

<<Audio`
<<Music`

(*Check your preception of dissonance as a function of frequency for two notes with different contents of harmonic partials.  The four cases assume both notes, one staying fixed and the other climbing in frequency, have the same harmonic structure. The first one is for two pure sine tones.  These take some time to calculate *)
amp1={2,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
amp2={2,3,0,0,0,0,0,0,0,0,0,0,0,0,0};
amp3={2,3,2,0,0,0,0,0,0,0,0,0,0,0,0};
amp4={1,3,2,1,.5,.3,.09,.5,.06,.3,.02,.07,.05,.03,.02};
f=200;
tmax=10;
fre[t_] =Apply[Plus,Table[amp1[[n]] Sin[2 Pi n f t],{n,1,15}]];
frechirp[t_] =Apply[Plus,Table[amp1[[n]] Sin[2 Pi n (f+f*t/tmax/2) t],{n,1,15}]];
Play[fre[t]+frechirp[t],{t,0,tmax},SampleRate->6*8192,SampleDepth->16] 
fre[t_] =Apply[Plus,Table[amp2[[n]] Sin[2 Pi n f t],{n,1,15}]];
frechirp[t_] =Apply[Plus,Table[amp2[[n]] Sin[2 Pi n (f+f*t/tmax/2) t],{n,1,15}]];
Play[fre[t]+frechirp[t],{t,0,tmax},SampleRate->6*8192,SampleDepth->16] 
fre[t_] =Apply[Plus,Table[amp3[[n]] Sin[2 Pi n f t],{n,1,15}]];
frechirp[t_] =Apply[Plus,Table[amp3[[n]] Sin[2 Pi n (f+f*t/tmax/2) t],{n,1,15}]];
Play[fre[t]+frechirp[t],{t,0,tmax},SampleRate->6*8192,SampleDepth->16] 
fre[t_] =Apply[Plus,Table[amp4[[n]] Sin[2 Pi n f t],{n,1,15}]];
frechirp[t_] =Apply[Plus,Table[amp4[[n]] Sin[2 Pi n (f+f*t/tmax/2) t],{n,1,15}]];
Play[fre[t]+frechirp[t],{t,0,tmax},SampleRate->6*8192,SampleDepth->16]

Dissonance Curve calculator
First, make a penalty function.  This should vanish for unison (no difference in frequency), rapidly rise to be irritating at 10-20 Hz, (beating),  and fall off with a rate that corresponds to distance on the basial membrane, falling to near zero when the distance becomes a critical bandwidth or more.

(* the penalty function; 0 at 0 Hz, and about 1/2 a critical bandwidth wide for 100 to 20000 Hz ceneter frequency (1/2 because it applies on either side of 0 frequency difference) *)
pen[x_,rng_]:=Sqrt[Abs[x]]Exp[-  (x/( 11+.02 rng))^2];
Plot[pen[x,100],{x,0,210},PlotStyle->Thick,PlotRange->All]

(*Predict dissonance as a function of frequency for different content of harmonic partials*)
(*Define a set of amplitudes, and apply the penalty to all frequency differences; sum up all penalties; plot as a function of f given f1.  Note many of the *)

f1=.;
(*Note the dips at many known musical intervals*)
base=200;
 amp4={1,3,2,1,1,1,.9,.5,.6,.3,.2,.2,.2,.1,.1};
dis[f1_,f_]:= Apply[Plus,Apply[Plus,Table[amp4[[i]]amp4[[j]]pen[i*f1-j*f,(i*f1+j*f)/2],{i,1,15},{j,1,15}]]];
a1=Plot[dis[base,f],{f,base,2*base+100},PlotRange->All,PlotStyle->Thick];
a2=Graphics[Table[{Line[{{base*i/j,5},{base*i/j,50}}],Text[Style[Rationalize[i/j],Large,Red],{base*i/j,60}]},{i,1,9},{j,Ceiling[i/2],i}]];
Show[{a1,a2}]